下载掌阅APP，畅读海量书库
立即打开
|
◆ |
第2章 金融理论模型 |
◆ |
本章用R语言深度解读了投资学理论和统计学理论在实际金融市场中的应用，包括4个基础理论模型，即资本资产定价模型、一元回归性线模型、多元回归线性模型、自回归模型，希望这些基础理论模型可以帮助读者，找到理解金融市场的方法。
问题
如何用R语言进行投资组合管理？
引言
伴随2016年中国金融交易市场的跌宕起伏，风险越来越不确定，利率持续走低，理财等无风险资产收益持续下降，唯有投资组合才能让我们的资产保值、增值。根据资本资产定价模型（CAPM），通过对金融数据的分析，构建投资组合，帮助我们在有效的市场中控制风险、稳定收益。
本节将深入浅出地介绍资本资产定价模型，从理论到建模，再到程序现实。资本资产定价模型反应的是资产的风险与期望收益之间的关系，风险越高，收益越高。当风险一样时，投资者会选择预期收益最高的资产；而预期收益一样时，投资者会选择风险最低的资产。 [1]
1952年，马科维茨（Markowitz）提出了投资组合选择理论，他认为最佳投资组合应当是风险厌恶特征的投资者的无差异曲线和资产的有效边界线的交点。投资者在选择资产时会在收益和风险之间做出平衡：当风险一样时，会选择预期收益最高的资产；而预期收益一样时，会选择风险最低的资产。
到1964年，威廉·夏普（William Sharp）、约翰·林特纳（John Lintner）与简·莫森（Jan Mossin）在马科维茨的基础上提出了单指数模型，将市场组合引入均值-方差模型，极大地简化了计算。他们认为获得了市场任意资组合的收益与某个共同因素之间有线性关系，最终将其发展为资本资产定价模型（Capital Asset Pricing Model，CAPM）。从马科维茨的投资组合选择理论，发展到资本资产定价模型经历了一个漫长的过程。
图2-1 投资组合选择示意图
简单一句话概括，资本资产定价模型的核心思想：资产价格取决于其获得的风险价格补偿。
资本资产定价模型，是基于一系列假设条件而成立的。但这些条件，可能并不符合现实的标准，资本资产定价模型也一度遭到质疑。
·资产可以无限分割。
·不存在交易成本和个人所得税。
·可以无限卖空。
·存在一种无风险利率，投资者在此利率水平下，可以无限制地贷出和借入任意数额的资金。
·投资者是价格接受者，市场是完全竞争的。
·投资者是理智的，通过比较资产的期望收益和方差来作出投资决策，在相同预期收益下会选择风险最小的资产。
·投资者在相同的投资期限出作出决策，而市场信息是公开免费的，并可以及时获得。
·投资者对市场中的经济变量有相同的预期，他们对任意资产的预期收益率、市场风险的看法是一致的。
资本资产定价模型的核心假设认为市场满足完全竞争、无摩擦和信息完全对称的条件，市场中的投资人都是Markowitz理论中的理性经济人。从真实的市场来看，这些假设条件基本都是不具备的。机构在用资本资产定价模型时候，也都是基于自己对市场的理解，在自己的优势端，尽量贴近假设条件，找出最优的投资组合方案，所以每个机构做的投资组合方案也都是不一样的。
由于涉及金融专业领域，有几个概念是我们应该提前知道的。
·风险资产：风险资产是指具有未来收益能力的资产，但收益率不确定且可能招致损失，比如股票、债券等。
·无风险资产：没有任何风险或者风险非常小的资产，有确定的收益率，并且不存在违约的风险。
·收益率：指从投资开始到投资结束时，所获得的投资回报率。
·无风险收益率：无风险资产，所产生的投资回报率。
·投资组合：由投资人或金融机构所持有的股票、债券、基金、衍生金融产品等组成的集合，目的在于分散风险。
·杠杆交易：就是利用小资金来进行数倍于原始金额的投资，以期望获取相对投资标的物波动的数倍收益率的盈利或亏损。
1.风险资产
对于风险资产来说，我们可以用预期收益和风险，通过二维的坐标来进行描述，如图2-2所示。
图2-2 风险资产
对图2-2的解释：
·x轴：风险；
·y轴：收益率；
·灰色区域：金融资产可投资区域；
·黑色线：有效投资边界；
·A点和B点：2个风险资产。
A和B有相同的x值，表示具有相同的风险。B点在A点上面，表示B的收益率高于A。对于理性的投资者来说，如果只在A点和B点之间做投资选择，那么大家会把所有钱都投资到B点，而不投资于A点。
2.无风险资产
在图2-3中，我们加入无风险资产，比较无风险资产和风险资产的关系。
对图2-3的解释：
·B点：1个风险资产，在有效投资边界上；
·C点：无风险资产，在y轴上；
·x轴：风险；
·y轴：收益率；
·灰色区域：金融资产为可投资区域；
·黑色线：有效投资边界。
C点为无风险资产，它的位置在图示的y轴上，这时x为0，即风险为0。我们可以把投资分配到C点或B点上。如果都投到C点，那么我们将获得的是r 0 部分的无风险收益；如果都投到B点，那么我们需要承担σ B 的风险，同时获得r B 的风险收益。如果我们把资金，一部分投资到B点对应的风险资产上，另一部分投资到C点对应的无风险资产上，那么将构成一个由B和C资产组成的投资组合，而且风险和收益部分将体现在B和C的连线上。
3.最优组合
那么，有没有最优的投资组合呢？收益最大、风险最小。下面就让我们来找出这个最优的组合M，如图2-4所示。
图2-3 无风险资产
图2-4 最优组合
对图2-4的解释：
·M点：最优组合的风险资产；
·B点：1个风险资产，在有效投资边界上；
·C点：无风险资产，在y轴上；
·x轴：风险；
·y轴：收益率；
·灰色区域：金融资产为可投资区域；
·黑色线：有效投资边界。
假设有最优的组合，在图2-4中M点处，当我们把C和M进行连线，使得CM的连线与灰色区域相切。从图上看，CM的连线会比任意的C点与可投资区域点的连线斜率都要大，比如C和B的连线。我们取CB的连线的延长线，在CB的延长线上找到，与M具有相同x的点B′，这时M与B′风险相同，M点在B′点的上面，所以M点的收益率大。也就是说，当风险相同的时候，我们都会选择收益率最大的资产。
不论从可投资区域中怎么选取，M点都是斜率最大的点，那么我们可以认为， M点为市场上各资产的最优的投资组合 。
对于最优的投资组合，其实不管投资者的收益风险的偏好是什么样子的，只要找到了最优的风险资产组合，再加上无风险的资产，就可以为投资者获得最佳的投资方案了。那么对于理性的投资者，如果发现了最优的组合，他们只会投资这个组合，这与收益和风险偏好无关。
M点构建的投资组合，一般是由所有可投资证券产品组成的，每种证券资产构成的比例，为证券的相对市值。无风险资产C，并没有包括在M中，人们都会选择CM的连接线进行投资，来构建最优的投资组合。
在实际的市场交易中，金融资产的价格会发生偏离，因为价格受市场的供需关系所影响，当价格发生偏离后，市场会自动修复会回均衡价格水平。
4.资本市场线
对于CM的连线，就是马科维茨提出的投资组合选择理论，即风险厌恶特征的投资者的无差异曲线和资产的有效边界线的交点。而CM这条线就叫资本市场线（Capital Market Line）。
资本市场线表明了有效组合的期望收益率和标准差之间的一种简单的线性关系。
资本市场线决定了证券的价格。因为资本市场线是证券有效组合条件下的风险与收益的均衡，如果脱离了这一均衡，则就会在资本市场线之外，形成另一种风险与收益的对应关系。
图2-5 资本市场线
5.投资组合构建
资本市场线，就是我们最优的投资组合，当我们发现这个投资组合，所有资金都会投到这个组合上。通过对无风险资产C和风险资产M分配不同的投资权重，我们可以配置出自己想要的风险和收益，同时可以利用金融工具来加杠杆放大风险和收益的范围，如图2-6所示。
我们把投资者分为风险厌恶型和风险激进型。
风险厌恶型对于资金安全有非常高的要求，不追求高收益但求本金安全，这些资金通常都是用来生活的。那么在为这些资金做资产配置方案的时候，可以把一部分资金配置到无风险资产上，同时将少量资金配置到M点的最优组合上，保证低风险并获得少量收益。
在图2-6中CM1点，如果配置50%的风险资产M和50%的无风险资产C，来实现投资组合，则公式如下：
CM 1 =0.5C+0.5M
对于风险激进型，他们对于资金有非常高的收益要求，本金可以部分或全部损失，这些资金通常都是“闲钱”，就是用来进行投资活动的。那么在为这些资金做资配置方案时，可以全部投到M上，再激进点，可以通过借钱、融资的方式，增加杠杆，把资金放大进行投资。这种操作风险会随着杠杆的放大剧增，当然同时你也会有更大的收益。
如图2-6中的CM 2 点，落在了CM的延长线上。我们可以配置150%的风险资产M，同时用50%的钱去抵押以无风险资产C的收益率去借钱。公式如下：
CM 2 =-0.5C+1.5M
6.风险和收益的关系
上面我们描述了风险和收益的关系，主要是从思路上定性介绍，没有进行定量描述，那么究竟风险和收益从数学上怎么进行定义呢？
图2-6 投资组合管理
图2-7 风险和收益关系
对图2-7的解释：
·M点：最优组合的风险资产；
·C点：无风险资产，在y轴上；
·r 0 ：无风险资产的收益率；
·r M ：M点的收益率；
·x轴：σ p 为风险资产的收益率的方差；
·y轴：r p 为收益率。
根据威廉·夏普所引入的均值-方差模型，极大地简化了计算，就是解决了公式计算的问题。用方差来刻画风险，建立收益和风险的一元线性关系。可以用下面公式来表示：
公式解释：
·E（r M ）：市场投资组合的预期收益率；
·r 0 ：无风险收益率；
·E（r M ）-r 0 ，市场投资组合的风险溢价；
·
：市场投资组合方差Var（r
M
）；
·A：风险厌恶水平。
有了公式，我们就明确知道风险和收益的定量关系，并且可以利用数据来进行计算。
对于市场的投资组合，风险溢价和市场投资组合的方差成线性关系。但对于单个资产来说，收益和风险是市场投资组合组成的一分部，受市场共同变化的影响。
1.单个资产风险溢价
对于单个资产的风险来说，在资本资产定价模型中，用β来进行表示。β是衡量单个金融资产与市场收益的共同变化程度，通过协方差来计算。单个资产的风险为，当前资产与投资组合收益率的协议差，除以投资组合收益率的方差。
单个资产的风险的计算公式：
单个资产的风险溢价的计算公式：
对公式的解释：
·E（r i ）：风险资产i的预期收益；
·E（r m ）：市场投资组合的预期收益；
·r f ：无风险资产收益；
·Cov（r i ，r m ）：风险资产收益率和市场投资组合收益率的协方差；
·Var（r m ）：市场投资组合的收益率的方差。
从公式可以看出，单个资产的风险溢价与市场投资组合M的风险溢价成正比，受β影响。
2.资本资产定价模型
资本资产定价模型，是现化金融学的基石理论。在上述假设条件下，可以推导出资本资产定价模型的具体公式。整个的推导过程，就是上面文章介绍的过程。从后人学习的角度看，这个理论是比较简单的，仅用到了简单地统计学知识，但是前人却花了很长的时间研究和探索。
由β值判断单个资产的风险，当β=1时，则说明当前资产与整个市场的趋势是完全保持一致的；当β为2时，代表高风险，其回报的变化将大于市场大盘的变化幅度；当β为0.5时，代表是低风险的资产配置。
3.2种风险
资本资产定价模型定义了2种风险，即系统性风险和非系统性风险。
系统性风险是由外部因素引起的风险，比如通货膨胀、GDP、重大政治事件等等。这一类事件对于资产收益率的影响不能通过组合本身来消除，所以这一类风险对于投资者来说是无法回避的。
非系统性风险是组合内部结构引起的风险，比如：A股与B股高度相关，A股的收益率出现大幅波动的时候，B股也会出现相似幅度的波动，从而发生波峰叠加或波谷叠加，就会增加整个组合的风险；反之，如果A与B为负相关，则A与B的波动就会相互抵消。这样，风险是由组合里的资产类型决定的，所以通过多样化分散的投资策略，无论在理论还是实际上，这种风险都是可以最小化甚至消除的。而在这个消除的过程中，整个投资组合的收益率是不会下降的。
4.2种收益
与风险相对应的是收益，我们承受了2种风险的同时，也获得了风险所带来的收益。一部分是与市场完全相关收益部分，即beta（β）收益；另一部分与市场不相关的收益部分，即alpha（α）收益。
beta收益相对容易获得，例如，你看好一个市场，以低成本买入对应市场的指数基金，等待市场上涨。alpha收益则比较难获得，alpha是体现投资水平的策略收益，是投资组合的实际期望收益与预期收益之间的差。计算alpha的公式为：
alpha是衡量投资人投资水平的，我们举个例子来说明。比如，市场收益率为14%，A证券的β=1.2，短期国债利率为6%，投资者对这只股票进行了交易，获得的实际收益为17%，那么我们怎么判断投资人的水平呢？
首先，求出A证券的预期收益率6%+1.2×（14-6）%=15.6%，再用投资者实际收益减去A证券预期收益17%-15.6%=1.4%。最后得到的1.4%就是alpha，表示以投资者能力，可以额外获得1.4%的收益。
5.资本资产定价模型的应用场景
进行组合投资分散风险：投资者可以按市场组合的构成比例分散持有多种风险资产，使持有的风险资产组合最大限度地接近市场组合，以达到消除非系统风险的目的。
调整收益风险比例：将无风险资产与风险资产市场组合进行再组合，以获得所希望的个性化的风险收益组合。
指数化投资：将资产配置在与某一指数相同的权重的投资方法，通过微调权重或成分，获得比指数更好的alpha。
资产定价：资本资产定价模型可以用来判断有价证券或其他金融资产的市场价格是否处于均衡水平、是否被高估或低估，以便通过套利活动获取超额收益。
基金购买：举一个贴近市场的例子，当我们要购买基金时，也可以用到资本资产定价模型，以帮助我们分析。比如，基金A的期望收益率12%，风险β=1，基金B期望收益率13%，β=1.5，市场期望收益率11%，无风险资产收益率r 0 =5%。那么哪只基金更值得买？
当你每天打开支付宝，看到里面的各种基金推荐，就会发现这是一个实际的问题。如果你学会了上述方法，按照资本资产定价模型的思路进行判断，其实购买基金就是求alpha，哪个基金的alpha高，就买哪个。
求alpha，我们就直接套用公式。
图2-8 基金alpha计算
基金A的alpha为1%，而基金B的alpha为-1%。结论就很明显，基金A的管理人能力很好，超额收益1%；而基金B的管理人，就差一些，盈利低于市场1%。所以，我们会投资基金A，而不会投资基金B。
上文用了大量的篇幅介绍了资本资产定价模型的原理，但它的程序实现其实是相当简单的。因为在R语言中，已经把资本资产定价模型相关的计算函数都封包好了，我们仅仅进行调用就能完成整个的计算过程。
R语言程序实现，我们主要会用到2个包，即quantmod包和PerformanceAnalytics包。对于R语言有哪些量化投资的包，可以参考1.1节。
·quantmod包用于下载数据；
·PerformanceAnalytics包用于进行各种评价指标计算。
我们设计一个应用场景，假如我们有10万美元要投资美国的股市，想获得比标普好（SP500）的投资收益，那么应该如何购买股票？
首先，要想清楚，我最终的目标是“比标普好的投资收益”。其次，我们基于资本资产定价模型理论基础，从投资组合角度思考投资策略，而不是技术指标的角度。比标普好，那么我们就需要以标普指数做为理想投资组合。然后，我们去市场上选择几个股票，分别计算出收益率、beta、alpha等指标，判断是否符合预期，反复测试，直到找到合适的股票或股票组合。
本节只是案例介绍，用于说明投资思路和方法，不构成任何的股票推荐。
本节所使用的系统环境：
·Win1064bit
·R：3.2.3 x86_64-w64-mingw32/x64 b4bit
从Yahoo下载IBM、GE（通用电器）、YHOO（Yahoo）3只股票，以及从2010年01月01日开始的日行情数据，同时下载标普指数（SP500）的日行情数据。
执行R语言程序：
# 加载程序包
> library(quantmod)
> library(PerformanceAnalytics)
# 从yahoo下载3只股票的数据，和SP500的数据
> getSymbols(c('IBM','GE','YHOO','^GSPC'), from = '2010-01-01')
# 打印前6行和后6行数据
> head(GSPC)
open high low close volume adjusted
2010-01-04 1116.56 1133.87 1116.56 1132.99 3991400000 1132.99
2010-01-05 1132.66 1136.63 1129.66 1136.52 2491020000 1136.52
2010-01-06 1135.71 1139.19 1133.95 1137.14 4972660000 1137.14
2010-01-07 1136.27 1142.46 1131.32 1141.69 5270680000 1141.69
2010-01-08 1140.52 1145.39 1136.22 1144.98 4389590000 1144.98
2010-01-11 1145.96 1149.74 1142.02 1146.98 4255780000 1146.98
> tail(GSPC)
open high low close volume adjusted
2016-12-20 2266.50 2272.56 2266.14 2270.76 3298780000 2270.76
2016-12-21 2270.54 2271.23 2265.15 2265.18 2852230000 2265.18
2016-12-22 2262.93 2263.18 2256.08 2260.96 2876320000 2260.96
2016-12-23 2260.25 2263.79 2258.84 2263.79 2020550000 2263.79
2016-12-27 2266.23 2273.82 2266.15 2268.88 1987080000 2268.88
2016-12-28 2270.23 2271.31 2249.11 2249.92 2392360000 2249.92
画出SP500的K线图，如图2-9所示。
> barChart(GSPC)
图2-9 SP500的K线图
对调整后的4个品种的价格进行合并。
> # 改列名
> names(IBM)<-c("open","high","low","close","volume","adjusted")
> names(GE)<-c("open","high","low","close","volume","adjusted")
> names(YHOO)<-c("open","high","low","close","volume","adjusted")
> names(GSPC)<-c("open","high","low","close","volume","adjusted")
# 数据合并
> dat=merge(IBM$adjusted,GE$adjusted,YHOO$adjusted,GSPC$adjusted)
> names(dat)<-c('IBM','GE','YHOO','SP500')
# 打印前6行
> head(dat)
IBM GE YHOO SP500
2010-01-04 112.2859 12.27367 17.10 1132.99
2010-01-05 110.9295 12.33722 17.23 1136.52
2010-01-06 110.2089 12.27367 17.17 1137.14
2010-01-07 109.8274 12.90920 16.70 1141.69
2010-01-08 110.9295 13.18724 16.70 1144.98
2010-01-11 109.7680 13.31435 16.74 1146.98
计算每日收益率，合并收益率到dat_ret。
> dat_ret=merge(IBM_ret,GE_ret,YHOO_ret,SP500_ret)
> names(dat_ret)<-c('IBM','GE','YHOO','SP500')
> head(dat_ret)
IBM GE YHOO SP500
2010-01-04 0.009681385 0.015111695 0.009445041 0.0147147759
2010-01-05 -0.012079963 0.005177994 0.007602339 0.0031156762
2010-01-06 -0.006496033 -0.005151320 -0.003482298 0.0005455205
2010-01-07 -0.003461515 0.051779935 -0.027373267 0.0040012012
2010-01-08 0.010034759 0.021538462 0.000000000 0.0028817272
2010-01-11 -0.010470080 0.009638554 0.002395150 0.0017467554
定义无风险收益率为4%，计算4个资产的平均年化收益率。
# 无风险收益率
> Rf<-0.04/12
# 计算平均年化收益率，平均年化标准差，平均年化Sharpe
> results<-table.AnnualizedReturns(dat_ret,Rf=Rf)
> results
IBM GE YHOO SP500
Annualized Return 0.0345 0.1108 0.1257 0.1055
Annualized Std Dev 0.1918 0.2180 0.3043 0.1555
Annualized Sharpe (Rf=84%) -2.8892 -2.3899 -1.6911 -3.3659
从平均年化收益率可以看到，这段时间标普的平均年化收益率有10%，已经是非常不错的成绩了，如果要想超过标普绝非易事。
我们再进行统计指标分析，从另外一个维度看一下数据的情况。
# 计算统计指标
> table.Stats(dat_ret)
IBM GE YHOO SP500
Observations 1760.0000 1760.0000 1760.0000 1760.0000
NAs 0.0000 0.0000 0.0000 0.0000
Minimum -0.0828 -0.0654 -0.0871 -0.0666
Quartile 1 -0.0060 -0.0065 -0.0098 -0.0039
Median 0.0002 0.0004 0.0005 0.0005
Arithmetic Mean 0.0002 0.0005 0.0007 0.0004
Geometric Mean 0.0001 0.0004 0.0005 0.0004
Quartile 3 0.0067 0.0077 0.0112 0.0053
Maximum 0.0567 0.1080 0.1034 0.0474
SE Mean 0.0003 0.0003 0.0005 0.0002
LCL Mean (0.95) -0.0004 -0.0001 -0.0002 0.0000
UCL Mean (0.95) 0.0008 0.0012 0.0015 0.0009
Variance 0.0001 0.0002 0.0004 0.0001
Stdev 0.0121 0.0137 0.0192 0.0098
Skewness -0.5876 0.3084 0.0959 -0.3514
Kurtosis 4.6634 4.7294 2.9990 4.0151
通过统计指标分析，每个资产有1760个样本点，没有NA值。在日最小收益率方面，YHOO为-0.0871。而在日最大收益率方面，GE为0.1080。算数平均、几何平均、方差、标准差都是YHOO最大。
画出IBM股票的日收益图（图2-10）、月收益图（图2-11），以及4个资产的累积收益率图（图2-12），并对4个资产做相关性分析。
图2-10 IBM股票每日收益
图2-11 IBM股票每月收益
图2-12 4个品种的累积收益率图
从图2-12中可以看出，GE和SP500的走势基本稳合，说明GE在从2010开始在跟着美国经济持续发展。YHOO从2013年年初到2015年年初大幅拉升，领先于SP500很多，说明这段时期YHOO所处的互联网行业，带来了非常大的市场红利；从2015年到2016年，又下跌很大，大起大落，受市场影响非常敏感。IBM大部分时间都处于SP500的下方，说明美国经济这几年的高速发展，并没有给IBM带来很大的发展空间。如果从我们的目标来说，“比标普好的投资收益”，那么我们只能选择GE或YHOO。
相关性分析
图2-13 相关性分析
对4个品种进行相关性分析，发现GE和SP500相关系数为0.78，是3只股票中相关度最高的。而YHOO是与其他3个品种走势最不一样的。
最后，以SP500为市场组合，分别计算出3只股票的alpha和beta。
# 计算alpha
> CAPM.alpha(dat_ret[,1:3],dat_ret[,4],Rf=Rf)
IBM GE YHOO
Alpha: SP500 -0.000752943 0.0003502332 0.0003944279
# 计算beta
> CAPM.beta(dat_ret[,1:3],dat_ret[,4],Rf=Rf)
IBM GE YHOO
Beta: SP500 0.8218135 1.098877 1.064844
在3只股票中，IBM的alpha是最小的，而且是负的，说明IBM落后于市场，买IBM不如直接买SP500更好。GE的Beta是最大的，在上升时期beta越大，获得的市场收益也会越大。从YHOO的Alpha和Beta上看，虽然与GE接近，但标准差、最大回撤等指标过大，会导致波动太大。
综上所述，配置部分GE和部分YHOO，就可以获得比标普好的收益，但由于GE和YHOO的beta都高于SP500，所以风险也会高于SP500，需要增加新的股票来分散风险，具体的定量分析将在后面的文章中再进行介绍。
最后，补充一些alpha和Beta的说明。alpha和beta的认知最早起源于一个股市的概念，是一个关于投资组合的收益率分解的问题。
·alpha：一般被认为是投资组合的超额收益，即管理人的能力。
·beta：市场风险，最初主要指股票市场的系统性风险。
alpha是平均实际回报和平均预期回报的差额。
·α>0：基金或股票的价格可能被低估，建议买入。
·α<0：基金或股票的价格可能被高估，建议卖空。
·α=0：基金或股票的价格准确反映其内在价值，未被高估也未被低估。
beta反映了单个证券与整体市场组合的联动性。
·β>1：攻击性，市场上升时涨幅大。
·β<1：防御性，市场下跌时跌幅小。
·β=1：中立性，与市场波动一致。
资本资产定价模型发展至今，已经有很长的时间了。金融理论在一直发展，继资本资产定价模型之后又一重要的理论突破是套利定价理论，我将在我的博客中进行介绍。
本节详细地介绍了资本资产定价模型的金融理论、推导过程，以及R语言实现。希望能给走在量化道路上的朋友带来入门的指引和帮助，也希望找到像我一样，通过IT转金融的人，让我们一起用IT技术+金融的思维在金融市场奋斗吧。
[1] 由于本节内容为非金融教材类，所以当出现与教课书不符的描述，请以教课书为准。本节力求用简化的语言，来介绍自资本资产定价模型的知识，同时配合R语言的实现。
问题
如何用R语言做回归分析？
引言
在我们的日常生活中，存在大量的具有相关性的事件，比如大气压和海拔高度，海拔越高大气压强越小；人的身高和体重，普遍来看越高的人体重也越重。还有一些可能存在相关性的事件，比如知识水平越高的人，收入水平越高；市场化的国家经济越好，则货币越强势，反之全球经济危机，黄金等避险资产则开始走强。
如果我们要研究这些事件，寻找不同变量之间的关系，就会用到回归分析。一元线性回归分析是处理两个变量之间线性相关关系最简单的模型。让我们一起发现生活中的规律吧。
回归分析（Regression Analysis）是用来确定2个或2个以上变量间关系的一种统计分析方法。如果在回归分析中，只包括一个自变量X和一个因变量Y，且它们的关系是线性的，那么这种回归分析称为一元线性回归分析。
回归分析属于统计学的基本模型，涉及统计学基础，需要掌握很多名词和知识点。
在回归分析中，变量有2类：因变量和自变量。因变量通常是指实际问题中所关心的指标，用Y表示：而自变量是影响因变量取值的一个变量，用X表示，如果有多个自变量则表示为X 1 ，X 2 ，…，X n 。
回归分析研究的主要步骤：
1）确定因变量Y与自变量X 1 ，X 2 ，…，X n 之间的定量关系表达式，即回归方程。
2）对回归方程的置信度进行检查。
3）判断自变量X n （n=1，2，…，m）对因变量的影响。
4）利用回归方程进行预测。
本节会根据回归分析的主要步骤（见图2-14），进行结构梳理，介绍一元线性回归模型的使用方法。
图2-14 回归分析的主要步骤
本节所使用的系统环境：
·Win1064bit；
·R：3.2.3 x86_64-w64-mingw32/x64 b4bit。
先让我们通过一个例子开始吧。下面用一组简单的数据来说明一元线性回归分析的数学模型的原理和公式。找出下面数据集中Y与X的定量关系。
数据集是2016年3月1日白天开盘的交易数据，为锌的2个期货合约的分钟线的价格数据，保存为文件zn.csv。数据集包括3列，其中Index列为时间，zn1.Close为ZN1604合约的1分钟线的报价数据，zn2.Close为ZN1605合约的1分钟线的报价数据。
打开zn.csv文件，数据集如下：
"Index","zn1.Close","zn2.Close" 2016-03-01 09:01:00,14075,14145 2016-03-01 09:02:00,14095,14160 2016-03-01 09:03:00,14095,14160 2016-03-01 09:04:00,14095,14165 2016-03-01 09:05:00,14120,14190 2016-03-01 09:06:00,14115,14180 2016-03-01 09:07:00,14110,14170 2016-03-01 09:08:00,14110,14175 2016-03-01 09:09:00,14105,14170 2016-03-01 09:10:00,14105,14170 2016-03-01 09:11:00,14120,14180 2016-03-01 09:12:00,14105,14170 2016-03-01 09:13:00,14105,14170 2016-03-01 09:14:00,14110,14175 2016-03-01 09:15:00,14105,14175 2016-03-01 09:16:00,14120,14185 2016-03-01 09:17:00,14125,14190 2016-03-01 09:18:00,14115,14185 2016-03-01 09:19:00,14135,14195 2016-03-01 09:20:00,14125,14190 2016-03-01 09:21:00,14135,14205 2016-03-01 09:22:00,14140,14210 2016-03-01 09:23:00,14140,14200 2016-03-01 09:24:00,14135,14205 2016-03-01 09:25:00,14140,14205 2016-03-01 09:26:00,14135,14205 2016-03-01 09:27:00,14130,14205 2016-03-01 09:28:00,14115,14180 # 省略
我们以zn1.Close列的价格为X，zn2.Close列的价格为Y，那么试试找到自变量X和因变量Y的关系的表达式。
为了直观起见，我们可以先画出一张散点图，以X为横坐标，Y为纵坐标，每个点对应一个X和一个Y。
# 加载类库
> library(zoo)
> library(xts)
# 设置环境变量
> options(stringsAsFactors = FALSE)
# 读取数据
> dat<-read.csv(f ile="zn.csv",sep=",")
# 转型为xts格式
> df<-xts(dat[,-1],order.by=as.POSIXct(dat[,1]))
# 数据集已存在df变量中
> head(df)
zn1.Close zn2.Close
2016-03-01 09:01:00 14075 14145
2016-03-01 09:02:00 14095 14160
2016-03-01 09:03:00 14095 14160
2016-03-01 09:04:00 14095 14165
2016-03-01 09:05:00 14120 14190
2016-03-01 09:06:00 14115 14180
# 分别给x,y赋值
> x<-as.numeric(df[,1])
> y<-as.numeric(df[,2])
# 画图
> plot(y~x+1)
图2-15 zn1.Close和zn2.Close散点图
从散点图上发现X和Y的排列基本是在一条直线附近，那么我们可以先假设X和Y的关系是线性，用公式表式为：
Y=a+b×X+c
·Y：因变量。
·X：自变量。
·a：截距。
·b：自变量系数。
·a+b×X：表示Y随X的变化而线性变化的部分。
·c：残差或随机误差，是其他一切不确定因素影响的总和，其值不可观测。假定c是符合均值为0，方差为σ 2 的正态分布，记作c~N（0，σ 2 ）。
对于上面的公式，称函数f（X）=a+b×X为一元线性回归函数，a为回归常数，b为回归系数，统称回归参数。X为回归自变量或回归因子，Y为回归因变量或响应变量。如果（X 1 ，Y 1 ），（X 2 ，Y 2 ）…（X n ，Y n ）是（X，Y）的一组观测值，则一元线性回归模型可表示为：
其中E（c i ）=0，var（c i ）=σ 2 ，i=1，2，...n
通过对一元线性回归模型的数学定义，接下来让我们利用数据集做回归模型的参数估计。
对于上面的公式，回归参数a和b是我们不知道的，我们需要用参数估计的方法来计算出a和b的值，而从得到数据集的X和Y的定量关系。我们的目标是要计算出一条直线，使直接线上每个点的Y值和实际数据的Y值之差的平方和最小，即（Y 1 实际-Y 1 预测） 2 +（Y 2 实际-Y 2 预测） 2 +……+（Y n 实际-Y n 预测） 2 的值最小。参数估计时，我们只考虑Y随X的线性变化的部分，而残差c是不可观测的，参数估计法并不需要考虑残差，对于残差的分析在后文中介绍。
令公式变形为a和b的函数Q（a，b），即（Y实际-Y测试）的平方和，变成（Y实际-（a+b×X））的平方和。
公式2-1回归参数变形公式：
通过最小二乘估计推导出a和b的求解公式。
公式2-2回归参数计算公式：
其中x和y的均值，计算方法如下：
有了这个公式，我们就可以求出a和b两个的回归参数的解了。
接下来，我们用R语言来实现对上面数据的回归模型的参数估计，R语言中可以用lm（）函数来实现一元线性回归的建模过程。
# 建立线性回归模型 > lm.ab<-lm(y ~ 1+x) # 打印参数估计的结果 > lm.ab Call: lm(formula = y ~ 1 + x) Coeff icients: (Intercept) x -349.493 1.029
如果你想动手来计算也可以自己实现公式。
# x均值 > Xm<-mean(x);Xm [1] 14034.82 # y均值 > Ym<-mean(y);Ym [1] 14096.76 # 计算回归系数 > b <- sum((x-Xm)*(y-Ym)) / sum((x-Xm)^2) ;b [1] 1.029315 # 计算回归常数 > a <- Ym - b * Xm;a [1] -349.493
回归参数a和b的计算结果与lm（）函数的计算结果是相同的。有了a和b的值，我们就可以画出这条近似的直线。
计算公式为：
Y= a + b * X = -349.493 + 1.029315 * X
画出回归线。
> plot(y~x+1) > abline(lm.ab)
图2-16 线性回归
这条直线是我们用数据拟合出来的，是一个近似的值。我们看到有些点在线上，有些点不在线上。那么要评价这条回归线拟合的好坏，我们就需要对回归模型进行显著性检验。
从公式2-2回归参数计算公式可知，在计算过程中并不一定要知道Y和X是否有线性相关的关系。如果不存相关关系，那么回归方程就没有任何意义了，如果Y和X是有相关关系的，即Y会随着X的变化而线性变化，这个时候一元线性回归方程才有意义。所以，我们需要用假设检验的方法，来验证相关性的有效性。
通常会采用3种显著性检验的方法。
·T检验法：T检验是检验模型某个自变量Xi对于Y的显著性，通常用P-value判断显著性，小于0.01时说明这个自变量Xi与Y相关关系显著。
·F检验法：F检验用于对所有的自变量X在整体上看对于Y的线性显著性，也是用P-value判断显著性，小于0.01时说明整体上自变量与Y相关关系显著。
·R 2 （R平方）相关系统检验法：用来判断回归方程的拟合程度，R 2 的取值在0，1之间，越接近1说明拟合程度越好。
在R语言中，上面列出的3种检验的方法都已被实现，我们只需要对结果进行解读。上文中，我们已经通过lm（）函数构建一元线性回归模型，然后可以summary（）函数来提取模型的计算结果。
> summary(lm.ab) # 计算结果
Call:
lm(formula = y ~ 1 + x)
Residuals:
Min 1Q Median 3Q Max
-11.9385 -2.2317 -0.1797 3.3546 10.2766
Coeff icients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.495e+02 7.173e+01 -4.872 2.09e-06 ***
x 1.029e+00 5.111e-03 201.390 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.232 on 223 degrees of freedom
Multiple R-squared: 0.9945, Adjusted R-squared: 0.9945
F-statistic: 4.056e+04 on 1 and 223 DF, p-value: < 2.2e-16
模型解读：
·Call：列出了回归模型的公式。
·Residuals：列出了残差的最小值点、1/4分位点、中位数点、3/4分位点，以及最大值点。
·Coefficients：表示参数估计的计算结果。
·Estimate：为参数估计列。Intercept行表示常数参数a的估计值，x行表示自变量x的参数b的估计值。
·Std.Error：为参数的标准差。
·t value：为t值，即T检验的值。
·Pr（>|t|）：表示P-value值，用于T检验判定，匹配显著性标记。
·显著性标记：***为非常显著，**为高度显著，*为显著，·为不太显著，没有记号为不显著。
·Residual standard error：表示残差的标准差，自由度为n-2。
·Multiple R-squared：为相关系数R 2 的检验，越接近1则越显著。
·Adjusted R-squared：为相关系数的修正系数，解决多元回归自变量越多，判定系数R 2 越大的问题。
·F-statistic：表示F统计量，自由度为（1，n-2），p-value用于F检验判定，匹配显著性标记。
通过查看模型的结果数据，我们可以发现通过T检验的截距和自变量x都是非常显著，通过F检验判断出整个模型的自变量是非常显著，同时R 2 的相关系数检验可以判断自变量和因变量是高度相关的。
最后，我们通过的回归参数的检验与回归方程的检验，得到一元线性回归方程为：
Y=-349.493+1.029315X
在得到的回归模型并进行显著性检验后，还要在做残差分析（预测值和实际值之间的差），检验模型的正确性，残差必须服从正态分布N（0，σ 2 ）。
我们可以自己计算数据残差，并进行正态分布检验。
# 残差
> y.res<-residuals(lm.ab)
# 打印前6条数据
> head(y.res)
1 2 3 4 5 6
6.8888680 1.3025744 1.3025744 6.3025744 5.5697074 0.7162808
# 正态分布检验
> shapiro.test(y.res)
Shapiro-Wilk normality test
data: y.res
W = 0.98987, p-value = 0.1164
# 画出残差散点图
> plot(y.res)
效果图如图2-17所示。
图2-17 残差散点图
对残差进行Shapiro-Wilk正态分布检验，W接近1，p-value>0.05，证明数据集符合正态分布。关于正态分布的介绍，可以参考《R的极客理想——高级开发篇》中1.4节的内容。
同时，我们也可以用R语言的工具生成4种用于模型诊断的图形，简化自己写代码计算的操作。
# 画图，回车展示下一张 > plot(lm.ab) Hit to see next plot: # 残差拟合图 Hit to see next plot: # 残差QQ图 Hit to see next plot: # 标准化的残差对拟合值 Hit to see next plot: # 标准化残差对杠杆值
图2-18 残差和拟合值对比图
对残差和拟合值作图，横坐标是拟合值，纵坐标是残差。残差和拟合值之间，数据点均匀分布在y=0两侧，呈现出随机的分布，红色线呈现出一条平稳的曲线并没有明显的形状特征，说明残差数据表现非常好。
残差QQ图，用来描述残差是否符合正态分布，如图2-19所示。图中的数据点按对角直线排列，趋于一条直线，并被对角直接穿过，直观上符合正态分布。对于近似服从正态分布的标准化残差，应该有95%的样本点落在[-2，2]区间内。
图2-19 残差QQ图
对标准化残差平方根和拟合值作图（见图2-20），横坐标是拟合值，纵坐标是标准化后的残差平方根。与残差和拟合值对比图（图2-18）的判断方法类似，数据随机分布，红色线呈现出一条平稳的曲线，无明显的形状特征。
图2-20 标准化残差平方根和拟合值对比图
对标准化残差和杠杆值作图（见图2-21），虚线表示的cooks距离等高线，通常用Cook距离度量的回归影响点。本图中没有出现红色的等高线，则说明数据中没有特别影响回归结果的异常点。
图2-21 标准残差和杠杆值对比图
如果想把4张图放在一起进行展示，可以改变画布布局，效果图如图2-22所示。
> par(mfrow=c(2,2)) > plot(lm.ab)
图2-22 合并4张图
在图2-22的4幅图中，每幅图上都有一些点被特别的标记出来了，这些点是可能存在的异常值点，如果要对模型进行优化，我们就可以从这些被标记出来的点入手。由于本次残差分析的结果是很好了，所以对于异常点的优化，可能并不能明显地提升模型的效果。
从图2-22中发现，索引编号为27和192的2个点在多幅图中出现。我们假设这2个点为异常点，从数据中去掉这2个点，再进行显著性检验和残差分析。
# 查看27和192
> df[c(27,192),]
zn1.Close zn2.Close
2016-03-01 09:27:00 14130 14205
2016-03-01 14:27:00 14035 14085
# 新建数据集，去掉27和192
> df2<-df[-c(27,192),]
回归建模和显著性检验。
> x2<-as.numeric(df2[,1])
> y2<-as.numeric(df2[,2])
> lm.ab2<-lm(y2 ~ 1+x2)
> summary(lm.ab2)
Call:
lm(formula = y2 ~ 1 + x2)
Residuals:
Min 1Q Median 3Q Max
-9.0356 -2.1542 -0.2727 3.3336 9.5879
Coeff icients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.293e+02 7.024e+01 -4.688 4.83e-06 ***
x2 1.028e+00 5.004e-03 205.391 < 2e-16 ***
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.117 on 221 degrees of freedom
Multiple R-squared: 0.9948, Adjusted R-squared: 0.9948
F-statistic: 4.219e+04 on 1 and 221 DF, p-value: < 2.2e-16
对比这次的显著性检验结果和之前结果，T检验，F检验和R 2 检验，并没有明显的效果提升，结果和我预想的是一样的。所以，通过残差分析和异常点分析，我认为模型是有效的。
最后，我们获得了一元线性回归方程的公式，就可以对数据进行预测了。比如，对给定X=x 0 时，计算出y 0 =a+bx 0 的值，并计算出置信度为1-α的预测区间。
当X=x 0 ，Y=y 0 时，置信度为1-α的预测区间为：
公式2-4预测区间：
公式2-5预测区间调整：
我们可以用R语言的predict（）函数来计算预测值y 0 ，和相应的预测区间。程序算法如下。
> newX<-data.frame(x=14040)
> lm.pred<-predict(lm.ab,newX,interval="prediction",level=0.95)
# 预测结果
> lm.pred
f it lwr upr
1 14102.09 14093.73 14110.44
当x 0 =14040时，在预测区间为0.95的概率时，y 0 的值为14102，预测区间为[14093.73，14110.44]。
我们通过图形来表示。
> plot(y~x+1)
> abline(lm.ab,col='red')
> points(rep(newX$x,3),y=lm.pred,pch=19,col=c('red','blue','green'))
图2-23 回归预测图
其中，红色点为y 0 的值，蓝色点为预测区间最小值，绿色点为预测区间最大值。
对于统计模型中最核心部分就在结果解读，本文介绍了一元回归模型的基本的建模过程和模型的详细解读方法。在我们掌握了这种方法以后，就可以更容易地理解和学习多元回归、非线性回归等更多的模型，并把这些模型应用到实际的工作中了。下一节将介绍用R语言解读多元线性回归模型，请继续阅读。
问题
如何用R语言做多元线性回归分析？
引言
在许多生活和工作的实际问题中，影响因变量的因素可能不止一个，比如对于知识水平越高的人，收入水平也越高，这样的一个论述。这其中可能包括了因为有更好的家庭条件，所以有了更好的教育；因为在一线城市发展，所以有了更好的工作机会；所处的行业赶上了大的经济上行周期等。这些共同的因素导致了较高的收入水平。这些规律是复杂的、多维度的，多元回归分析方法更适合解读生活的规律。
对比一元线性回归，多元线性回归是用来确定2个或2个以上变量间关系的统计分析方法。多元线性回归的基本的分析方法与一元线性回归方法是类似的，我们首先需要选取多元数据集并定义数学模型，然后进行参数估计，对估计出来的参数进行显著性检验、残差分析、异常点检测，最后确定回归方程并进行模型预测。
由于多元回归方程有多个自变量，区别于一元回归方程，有一项很重要的操作就是自变量的优化，挑选出相关性最显著的自变量，同时去除不显著的自变量。在R语言中，有很方便地用于优化自变量的函数，可以很好地帮助我们改进回归模型。下面就开始多元线性回归的建模过程。
本节所使用的系统环境：
·Win1064bit
·R：3.2.3 x86_64-w64-mingw32/x64 b4bit
做过商品期货研究的人，都知道黑色系品种是具有产业链上下游的关系。铁矿石是炼钢的原材料，焦煤和焦炭是炼钢的能源资源，热卷即热轧卷板是以板坯为原料经加热后制成的钢板，螺纹钢是表面带肋的钢筋。
由于产业链的关系，假设我们想要预测螺纹钢的价格，那么影响螺纹钢价格的因素可能会涉及原材料、能源资源和同类材料等。比如，我们直观上认为，如果铁矿石价格上涨，螺纹钢就应该要涨价了。那么当铁矿石1吨涨1000元后，螺纹钢1吨涨多少钱？靠经验就可能估计不准，需要定量地做计算。
1.数据集和数学模型
先从数据开始介绍，这次的数据集，我选择的是黑色系的商品期货，包括了大连期货交易所的焦煤（JM）、焦炭（J）、铁矿石（I），上海期货交易所的螺纹钢（RU）和热卷（HC）。
数据集为2016年3月15日当日白天开盘的交易数据，为黑色系的5个期货合约的分钟线的价格数据，保存为文件future.csv。
# 加载类库
> library(xts)
> library(reshape2)
> library(ggplot2)
# 设置环境变量
> options(stringsAsFactors = FALSE)
# 读取数据
> dat<-read.csv(f ile="future.csv",sep=",")
# 转型为xts格式
> df<-xts(dat[,-1],order.by=as.Date(dat[,1]))
# 数据集已存在df变量中
> head(df,20)
x1 x2 x3 x4 y
2016-03-15 09:01:00 754.5 616.5 426.5 2215 2055
2016-03-15 09:02:00 752.5 614.5 423.5 2206 2048
2016-03-15 09:03:00 753.0 614.0 423.0 2199 2044
2016-03-15 09:04:00 752.5 613.0 422.5 2197 2040
2016-03-15 09:05:00 753.0 615.5 424.0 2198 2043
2016-03-15 09:06:00 752.5 614.5 422.0 2195 2040
2016-03-15 09:07:00 752.0 614.0 421.5 2193 2036
2016-03-15 09:08:00 753.0 615.0 422.5 2197 2043
2016-03-15 09:09:00 754.0 615.5 422.5 2197 2041
2016-03-15 09:10:00 754.5 615.5 423.0 2200 2044
2016-03-15 09:11:00 757.0 616.5 423.0 2201 2045
2016-03-15 09:12:00 756.0 615.5 423.0 2200 2044
2016-03-15 09:13:00 755.5 615.0 423.0 2197 2042
2016-03-15 09:14:00 755.5 615.0 423.0 2196 2042
2016-03-15 09:15:00 756.0 616.0 423.5 2200 2045
2016-03-15 09:16:00 757.5 616.0 424.0 2205 2052
2016-03-15 09:17:00 758.5 618.0 424.0 2204 2051
2016-03-15 09:18:00 759.5 618.5 424.0 2205 2053
2016-03-15 09:19:00 759.5 617.5 424.5 2206 2053
2016-03-15 09:20:00 758.5 617.5 423.5 2201 2050
数据集包括6列：
·索引：时间；
·x1：焦炭（j1605）合约的1分钟线的报价数据；
·x2：焦煤（jm1605）合约的1分钟线的报价数据；
·x3：铁矿石（i1605）合约的1分钟线的报价数据；
·x4：热卷（hc1605）合约的1分钟线的报价数据；
·y：螺纹钢（rb1605）合约的1分钟线的报价数据。
假设螺纹钢的价格与其他4个商品的价格有线性关系，那么我们建立以螺纹钢为因变量，以焦煤、焦炭、铁矿石和热卷的为自变量的多元线性回归模型。用公式表示为：
·y：因变量，螺纹钢；
·x1：自变量，焦煤；
·x2：自变量，焦炭；
·x3：自变量，铁矿石；
·x4：自变量，热卷；
·a：截距；
·b、c、d、e：自变量系数；
·ε：残差，是其他一切不确定因素影响的总和，其值不可观测。假定ε服从正态分布N（0，σ^2）。
通过对多元线性回归模型的数学定义，接下来让我们利用数据集做多元回归模型的参数估计。
2.回归参数估计
上面公式中，回归参数a、b、c、d、e都是我们不知道的，参数估计就是通过数据来估计出这些参数，从而确定自变量和因变量之前的关系。我们的目标是要计算出一条直线，使直线上每个点的Y值和实际数据的Y值之差的平方和最小，即（Y1实际-Y1预测）^2+（Y2实际-Y2预测）^2+……+（Yn实际-Yn预测）^2的值最小。参数估计时，我们只考虑Y随X自变量的线性变化的部分，而残差ε是不可观测的，参数估计法并不需要考虑残差。
类似于一元线性回归，我们用R语言来实现对数据的回归模型的参数估计，用lm（）函数来实现多元线性回归的建模过程。
# 建立多元线性回归模型 > lm1<-lm(y~x1+x2+x3+x4,data=df) # 打印参数估计的结果 > lm1 Call: lm(formula = y ~ x1 + x2 + x3 + x4, data = df) Coeff icients: (Intercept) x1 x2 x3 x4 212.8780 0.8542 0.6672 -0.6674 0.4821
这样我们就得到了y和x关系的方程。
3.回归方程的显著性检验
参考一元线性回归的显著性检验，多元线性回归的显著性检验，同样是需要经过T检验，F检验，以及R^2（R平方）相关系统检验。在R语言中这三种检验的方法都已被实现，我们只需要把结果解读，我们可以summary（）函数来提取模型的计算结果。
> summary(lm1)
Call:
lm(formula = y ~ x1 + x2 + x3 + x4, data = df)
Residuals:
Min 1Q Median 3Q Max
-4.9648 -1.3241 -0.0319 1.2403 5.4194
Coeff icients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 212.87796 58.26788 3.653 0.000323 ***
x1 0.85423 0.10958 7.795 2.50e-13 ***
x2 0.66724 0.12938 5.157 5.57e-07 ***
x3 -0.66741 0.15421 -4.328 2.28e-05 ***
x4 0.48214 0.01959 24.609 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.028 on 221 degrees of freedom
Multiple R-squared: 0.9725, Adjusted R-squared: 0.972
F-statistic: 1956 on 4 and 221 DF, p-value: < 2.2e-16
·T检验：所有自变量都是非常显著***；
·F检验：同样是非常显著，p-value<2.2e-16；
·调整后的R 2 ：相关性非常强为0.972。
最后，我们通过的回归参数的检验与回归方程的检验，得到最后多元线性回归方程：
即
4.残差分析和异常点检测
在得到的回归模型进行显著性检验后，还要再做残差分析（预测值和实际值之间的差），检验模型的正确性，残差必须服从正态分布N（0，σ^2）。直接用plot（）函数生成4种用于模型诊断的图形，进行直观地分析。
> par(mfrow=c(2,2)) > plot(lm1)
·残差和拟合值（左上），残差和拟合值之间数据点均匀分布在y=0两侧，呈现出随机的分布，红色线呈现出一条平稳的曲线并没有明显的形状特征。
·残差QQ图（右上），数据点按对角直线排列，趋于一条直线，并被对角直接穿过，直观上符合正态分布。
·标准化残差平方根和拟合值（左下），数据点均匀分布在y=0两侧，呈现出随机的分布，红色线呈现出一条平稳的曲线并没有明显的形状特征。
·标准化残差和杠杆值（右下），没有出现红色的等高线，则说明数据中没有特别影响回归结果的异常点。
结论：没有明显的异常点，残差符合假设条件。
图2-24 多元回归残差检验
5.模型预测
我们得到了多元线性回归方程的公式，就可以对数据进行预测了。我们可以用R语言的predict（）函数来计算预测值y 0 和相应的预测区间，并把实际值和预测值一起可视展示。
> par(mfrow=c(1,1)) #设置画面布局
# 预测计算
> dfp<-predict(lm1,interval="prediction")
# 打印预测值
> head(dfp,10)
f it lwr upr
2014-03-21 3160.526 3046.425 3274.626
2014-03-24 3193.253 3078.868 3307.637
2014-03-25 3240.389 3126.171 3354.607
2014-03-26 3228.565 3114.420 3342.710
2014-03-27 3222.528 3108.342 3336.713
2014-03-28 3262.399 3148.132 3376.666
2014-03-31 3291.996 3177.648 3406.344
2014-04-01 3305.870 3191.447 3420.294
2014-04-02 3275.370 3161.018 3389.723
2014-04-03 3297.358 3182.960 3411.755
# 合并数据
> mdf<-merge(df$y,dfp)
# 画图，如图2-25所示
> draw(mdf)
图2-25 模型预测分析
图例说明
·y：实际价格，红色线；
·fit：预测价格，绿色线；
·lwr：预测最低价，蓝色线；
·upr：预测最高价，紫色线。
从图中看出，实际价格y和预测价格fit，在大多数的时候都是很贴近的。我们的一个模型就训练好了！
上文中，我们已经很顺利地发现了一个非常不错的模型。如果要进行模型优化，可以用R语言中update（）函数进行模型的调整。我们首先检查一下每个自变量x1、x2、x3、x4和因变量y之间的关系，其配对展示图如图2-26所示。
# 各变量之间的关系 > pairs(as.data.frame(df))
从图2-26中我们可以发现，x2与y的关系可能是最偏离线性的。那么，我们尝试对多元线性回归模型进行调整，从原模型中去掉x2变量。
# 模型调整
> lm2<-update(lm1, .~. -x2)
> summary(lm2)
Call:
lm(formula = y ~ x1 + x3 + x4, data = df)
Residuals:
Min 1Q Median 3Q Max
-6.0039 -1.3842 0.0177 1.3513 4.8028
Coeff icients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 462.47104 34.26636 13.50 < 2e-16 ***
x1 1.08728 0.10543 10.31 < 2e-16 ***
x3 -0.40788 0.15394 -2.65 0.00864 **
x4 0.42582 0.01718 24.79 < 2e-16 ***
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.142 on 222 degrees of freedom
Multiple R-squared: 0.9692, Adjusted R-squared: 0.9688
F-statistic: 2330 on 3 and 222 DF, p-value: < 2.2e-16
图2-26 配对展示
当把自变量x2去掉后，自变量x3的T检验值反而变大了，同时Adjusted R-squared变小了，所以我们这次调整是有问题的。
通过生产和原材料的内在逻辑分析，焦煤与焦炭属于上下游关系。焦煤是生产焦炭的一种原材料，焦炭是焦煤与其他炼焦煤经过配煤焦化形成的产品，一般生产1吨焦炭需要1.33吨炼焦煤，其中焦煤至少占30%。
我们把焦煤和焦炭的关系改变一下，增加x1*x2的关系并匹配到模型，看看效果。
# 模型调整
> lm3<-update(lm1, .~. + x1*x2)
> summary(lm3)
Call:
lm(formula = y ~ x1 + x2 + x3 + x4 + x1:x2, data = df)
Residuals:
Min 1Q Median 3Q Max
-4.8110 -1.3501 -0.0595 1.2019 5.3884
Coeff icients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7160.32231 7814.50048 0.916 0.361
x1 -8.45530 10.47167 -0.807 0.420
x2 -10.58406 12.65579 -0.836 0.404
x3 -0.64344 0.15662 -4.108 5.63e-05 ***
x4 0.48363 0.01967 24.584 < 2e-16 ***
x1:x2 0.01505 0.01693 0.889 0.375
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.029 on 220 degrees of freedom
Multiple R-squared: 0.9726, Adjusted R-squared: 0.972
F-statistic: 1563 on 5 and 220 DF, p-value: < 2.2e-16
从结果中发现，增加了x1*x2列后，原来的x1、x2和Intercept的T检验都不显著。继续调整模型，从模型中去掉x1，x2两个自变量。
# 模型调整
> lm4<-update(lm3, .~. -x1-x2)
> summary(lm4)
Call:
lm(formula = y ~ x3 + x4 + x1:x2, data = df)
Residuals:
Min 1Q Median 3Q Max
-4.9027 -1.2516 -0.0167 1.2748 5.8683
Coeff icients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.950e+02 1.609e+01 43.183 < 2e-16 ***
x3 -6.284e-01 1.530e-01 -4.108 5.61e-05 ***
x4 4.959e-01 1.785e-02 27.783 < 2e-16 ***
x1:x2 1.133e-03 9.524e-05 11.897 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.035 on 222 degrees of freedom
Multiple R-squared: 0.9722, Adjusted R-squared: 0.9718
F-statistic: 2588 on 3 and 222 DF, p-value: < 2.2e-16
从调整后的结果来看，效果还不错。不过，也并没有比最初的模型有所提高。
对于模型调整的过程，在我们手动调整测试时，一般都会基于业务知识来操作。如果是按照数据指标来计算，我们可以用R语言中提供的逐步回归的优化方法，通过AIC指标来判断是否需要参数优化。
#对lm1模型做逐步回归
> step(lm1)
Start: AIC=324.51
y ~ x1 + x2 + x3 + x4
Df Sum of Sq RSS AIC
908.8 324.51
- x3 1 77.03 985.9 340.90
- x2 1 109.37 1018.2 348.19
- x1 1 249.90 1158.8 377.41
- x4 1 2490.56 3399.4 620.65
Call:
lm(formula = y ~ x1 + x2 + x3 + x4, data = df)
Coeff icients:
(Intercept) x1 x2 x3 x4
212.8780 0.8542 0.6672 -0.6674 0.4821
通过计算AIC指标，lm1的模型AIC最小时为324.51，每次去掉一个自变量都会让AIC的值变大，所以我们还是不调整比较好。
对刚才的lm3模型做逐步回归的模型调整。
#对lm3模型做逐步回归
> step(lm3)
Start: AIC=325.7 #当前AIC
y ~ x1 + x2 + x3 + x4 + x1:x2
Df Sum of Sq RSS AIC
- x1:x2 1 3.25 908.8 324.51
905.6 325.70
- x3 1 69.47 975.1 340.41
- x4 1 2487.86 3393.5 622.25
Step: AIC=324.51 #去掉x1*x2项的AIC
y ~ x1 + x2 + x3 + x4
Df Sum of Sq RSS AIC
908.8 324.51
- x3 1 77.03 985.9 340.90
- x2 1 109.37 1018.2 348.19
- x1 1 249.90 1158.8 377.41
- x4 1 2490.56 3399.4 620.65
Call:
lm(formula = y ~ x1 + x2 + x3 + x4, data = df)
Coeff icients:
(Intercept) x1 x2 x3 x4
212.8780 0.8542 0.6672 -0.6674 0.4821
通过AIC的判断，去掉x1*x2项后AIC最小，最后的检验结果告诉我们，还是原初的模型是最好的。
最后，我们用上面5个期货合约的日K线数据测试一下，找到多元回归关系。加载黑色系5个品种的期货数据，数据分别存储在文件中，期货品种与文件名的对应关系如表2-1所示。
表2-1 期货品种与文件名对应关系
接下来，用R语言进行读取数据和数据格式转换。
# 定义读文件函数
> dailyData <- function(f ile) {
+ df <- read.table(f ile = f ile, header = FALSE,sep = ',', na.strings = 'NULL')
+ names(df)<-c('date','price')
+ return(df)
+ }
# 转型为xts类型
> toXts<-function(data,format='%Y-%m-%d %H:%M:%S'){
+ df<-subset(data,select=-c(date))
+ xts(df,order.by=strptime(data$date, format))
+ }
# 分别从文件中读取数据
> x1<-toXts(dailyData(f ile='j_daily.csv'),'%Y%m%d')
> x2<-toXts(dailyData(f ile='jm_daily.csv'),'%Y%m%d')
> x3<-toXts(dailyData(f ile='i_daily.csv'),'%Y%m%d')
> x4<-toXts(dailyData(f ile='hc_daily.csv'),'%Y%m%d')
> y<-toXts(dailyData(f ile='rb_daily.csv'),'%Y%m%d')
# 合并数据集
> df<-na.omit(merge(x1,x2,x3,x4,y))
# 对列进行重命名
> names(df)<-c('x1','x2','x3','x4','y')
对上面5个品种的日K线数据进行多元回归分析建模。
# 多元回归分析建模
> lm9<-lm(y~x1+x2+x3+x4,data=df) # 日K线数据
> summary(lm9)
Call:
lm(formula = y ~ x1 + x2 + x3 + x4, data = df)
Residuals:
Min 1Q Median 3Q Max
-173.338 -37.470 3.465 32.158 178.982
Coeff icients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 386.33482 31.07729 12.431 < 2e-16 ***
x1 0.75871 0.07554 10.045 < 2e-16 ***
x2 -0.62907 0.14715 -4.275 2.24e-05 ***
x3 1.16070 0.05224 22.219 < 2e-16 ***
x4 0.46461 0.02168 21.427 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 57.78 on 565 degrees of freedom
Multiple R-squared: 0.9844, Adjusted R-squared: 0.9843
F-statistic: 8906 on 4 and 565 DF, p-value: < 2.2e-16
查看数据集的基本统计信息。
> summary(df)
Index x1 x2
Min. :2014-03-21 00:00:00 Min. : 606.5 Min. :494.0
1st Qu.:2014-10-21 06:00:00 1st Qu.: 803.5 1st Qu.:613.1
Median :2015-05-20 12:00:00 Median : 939.0 Median :705.8
Mean :2015-05-21 08:02:31 Mean : 936.1 Mean :695.3
3rd Qu.:2015-12-16 18:00:00 3rd Qu.:1075.0 3rd Qu.:773.0
Max. :2016-07-25 00:00:00 Max. :1280.0 Max. :898.0
x3 x4 y
Min. :284.0 Min. :1691 Min. :1626
1st Qu.:374.1 1st Qu.:2084 1st Qu.:2012
Median :434.0 Median :2503 Median :2378
Mean :476.5 Mean :2545 Mean :2395
3rd Qu.:545.8 3rd Qu.:2916 3rd Qu.:2592
Max. :825.0 Max. :3480 Max. :3414
最后，生成期货日K线的5个品种的配对关系展示图，如图2-27所示。
> pairs(as.data.frame(df))
图2-27 黑色商品期货的相关关系
对于黑色系的5个品种的日K线数据，同样具有非常强的相关关系，那么我们就可以把这个结论应用到实际的交易中了。
本文通过多元回归的统计分析方法，介绍多元回归在金融市场的基本应用。我们通过建立因变量和多个自变量的模型，从而发现生活中更复杂的规律，并建立有效的验证指标。让我用技术的优势，去更多的发现金融市场的规律吧。
问题
如何用R语言进行时间序列分析？
引言
时间序列是金融分析中常用到的一种数据格式，自回归模型是分析时间序列数据的一种基本的方法。通过建立自回归模型，找到数据自身周期性的规律，从而帮助我们理解金融市场的变化。
在时间序列分析中，有很多常用的模型，包括AR、MA、ARMA、ARIMA、ARCH、GARCH，这些模型的主要区别是适用条件不同，且是层层递进的，后面的一个模型解决了前一个模型的某个固有问题。本节从AR模型开始，将对时间序列分析体系进行完整的介绍，并用R语言进行模型实现。
自回归模型（Autoregressive Model，AR模型），是统计上一种处理时间序列的方法，用来描述当前值与历史值之间的关系，用变量自身的历史时间数据对自身进行预测，自回归模型必须满足平稳性的要求。比如，时间序列数据集X的历史各期数据从X 1 至X t-1 ，假设它们为线性关系，可以对当期Xt的表现进行预测。X的当期值等于一个或数个落后期的线性组合，加常数项，加随机误差。
p阶自回归过程的公式定义：
字段解释：
·X t ：当期的X的表现；
·c：常数项；
·p：阶数，i为从1到p的值；
·φ i ：自相关系数；
·t：时间周期；
·ε t ：均值为0，标准差为δ的随机误差，同时δ是独立于t的。
对于一阶自回模型，用AR（1）来表示，简化后的公式为：
自回归是从线性回归分析中发展而来，只是把自变量x对因变量y的分析，变成自变量x对自身的分析。如果你需要了解线性回归的知识，请参考2.2节。
1.自回归模型的限制
自回归模型是用自身的数据来进行预测，但是这种方法受到一定的限制：
·必须具有平稳性，平稳性要求随机过程的随机特征不随时间变化。
·必须具有自相关性，如果自相关系数（φ i ）小于0.5，则不宜采用，否则预测结果极不准确。
·自回归只适用于预测与自身前期相关的现象，即受自身历史因素影响较大的现象。对于受其他因素影响的现象，不宜采用自回归，可以改用向量自回归模型。
2.平稳性时间序列的特点
平稳性要求产生时间序列Y的随机过程的随机特征不随时间变化，则称过程是平稳的；假如该随机过程的随机特征随时间变化，则称过程是非平稳的。
平稳性是由样本时间序列所得到的拟合曲线，在未来的一段期间内能顺着现有的形态能一直地延续下去；如果数据非平稳，则说明样本拟合曲线的形态不具有延续的特点，也就是说拟合出来的曲线将不符合当前曲线的形态。
·随机变量Y t 的均值和方差均与时间t无关；
·随机变量Y t 和Y s 的协方差只与时间差（步长）t-s有关；
·对于平稳时间序列在数学上有比较丰富的处理手段，非平稳的时间序列通过差分等手段转化为平稳时间序列处理。
了解了自回归模型的定义，我们就可以用R语言来模拟一下自回归模型的构建和计算过程。
本文所使用的系统环境：
·Win1064bit
·R：3.2.3 x86_64-w64-mingw32/x64 b4bit
首先，我们生成一个随机游走的数据集，满足平稳性的要求。
# 随机游走的数据集
> set.seed(0)
> x<-w<-rnorm(1000) # 生成符合正态分布N(0,1)的数据
> for(t in 2:1000) x[t]<-x[t-1]+w[t]
> tsx<-ts(x) # 生成ts时间序列的数据集
# 查看数据集
> head(tsx,15)
[1] 1.2629543 0.9367209 2.2665202 3.5389495 3.9535909 2.4136409
[7] 1.4850739 1.1903534 1.1845862 3.5892396 4.3528331 3.5538238
[13] 2.4061668 2.1167053 1.8174901
> plot(tsx) # 生成可视化图形
> a<-ar(tsx);a # 进行自回归建模
Call:
ar(x = tsx)
Coeff icients:
1
0.9879
Order selected 1 sigma^2 estimated as 1.168
图2-28 随机游走的数据集
自相关系数为0.9879，这是一个非常强的自相关性，所以上述的数列符合自相关的特性。
R语言中ar（）函数提供了多种自相关系数的估计，包括“yule-walker”“burg”“ols”“mle”“yw”，默认是用yule-walker方法，常用的方法还有最小二乘法（ols），极大似然法（mle）。
我们用最小二乘法来进行参数估计。
> b<-ar(tsx,method = "ols");b
Call:
ar(x = tsx, method = "ols")
Coeff icients:
1
0.9911
Intercept: -0.017 (0.03149)
Order selected 1 sigma^2 estimated as 0.9906
用最小二乘法的计算结果，则自相关系统数为0.9911，截距为-0.017。只有使用最小二乘法进行参数估计的时候，才会有截距。
我们用极大似然法来进行参数估计。
> d<-ar(tsx,method = "mle");d
Call:
ar(x = tsx, method = "mle")
Coeff icients:
1
0.9904
Order selected 1 sigma^2 estimated as 0.9902
用极大似然法计算结果，则自相关系统数为0.9904。对于上面3种估计方法，自相关系数的值都是很接近的。
在上面的例子中，我们默认是用一阶的自回归模型AR（1），进行程序实现的。在实际应用中，自回归模型AR时间序列的阶数P是未知的，必须根据实际数据来决定，就要对AR模型定阶数。常的方法就是利用自相关函数（ACF）和偏自相关函数（PACF）来确定自回归模型的阶数。在ACF/PACF不能确定的情况下，还需要用AIC（Aikaike info Criterion）、BIC（Bayesian information criterion）的信息准则函数来确定阶数。
自回归模型的建模过程，是通过确定阶数、参数估计、再次确定阶数的方法进行判断。自相关函数ACF可确定采用自回归模型是否合适。如果自相关函数具有拖尾性，则AR模型为合适模型。偏自相关函数PACF用来确定模型的阶数，如果从某个阶数之后，偏自相关函数的值都很接近0，则取相应的阶数作为模型阶数，偏自相关函数通过截尾性确定阶数。
1.自相关函数ACF（autocorrelation function）
将一个有序的随机变量序列与其自身相比较，这就是自相关函数在统计学中的定义。每个不存在相位差的序列都与其自身相似，即在此情况下，自相关函数值最大。如果序列中的组成部分相互之间存在相关性（不再是随机的），则由以下相关值方程所计算的值不再为零，这样的组成部分为自相关。
自相关函数反映了同一序列在不同时序的取值之间的相关程序。
ACF的公式为：
字段解释
·P k ，为ACF的标准误差；
·t，为数据集的长度；
·k，为滞后，取值从1到t-1，表示相距k个时间间隔的序列值之间的相关性；
·Y t ，为样本在t时期的值；
·Y t-k ，为样本在t-k时期的值；
·μ，为样本的均值。
所得的自相关值P k 的取值范围为[-1，1]，1为最大正相关值，-1则为最大负相关值，0为不相关。
根据上面公式，我们也可以手动计算出tsx数据集的ACF值。
> u<-mean(tsx) #均值 > v<-var(tsx) #方差 > # 1阶滞后 > p1<-sum((x[1:length(tsx)-1]-u)*(x[2:length(tsx)]-u))/((length(tsx)-1)*v);p1 [1] 0.9878619 > # 2阶滞后 > p2<-sum((x[1:(length(tsx)-2)]-u)*(x[3:length(tsx)]-u))/((length(tsx)-1)*v);p2 [1] 0.9760271 > # 3阶滞后 > p3<-sum((x[1:(length(tsx)-3)]-u)*(x[4:length(tsx)]-u))/((length(tsx)-1)*v);p3 [1] 0.9635961
同时，我们可以用R语言中的acf（）函数来计算，打印前30个滞后的ACF值。
> acf(tsx)$acf
, , 1
[,1]
[1,] 1.0000000
[2,] 0.9878619
[3,] 0.9760271
[4,] 0.9635961
[5,] 0.9503371
[6,] 0.9384022
[7,] 0.9263075
[8,] 0.9142540
[9,] 0.9024862
[10,] 0.8914740
[11,] 0.8809663
[12,] 0.8711005
[13,] 0.8628609
[14,] 0.8544984
[15,] 0.8462270
[16,] 0.8384758
[17,] 0.8301834
[18,] 0.8229206
[19,] 0.8161523
[20,] 0.8081941
[21,] 0.8009467
[22,] 0.7942255
[23,] 0.7886249
[24,] 0.7838154
[25,] 0.7789733
[26,] 0.7749697
[27,] 0.7709313
[28,] 0.7662547
[29,] 0.7623381
[30,] 0.7604101
[31,] 0.7577333
比较前3个值的计算结果，与我们自己的计算结果是一样的，同时可以用R语言进行可视化输出。
> acf(tsx)
从图2-29中看出，数据的ACF为拖尾，存在很严重的自相关性。接下来，这时候我们用偏自相关函数确定一下AR的阶数。
图2-29 ACF检验
2.偏自相关函数（PACF）（partial autocorrelation function）
偏自相关函数是由自相关函数推到而来。对于一个平稳AR（p）模型，求出滞后k自相关系数p（k）时，实际上得到并不是x（t）与x（t-k）之间单纯的相关关系。因为x（t）同时还会受到中间k-1个随机变量x（t-1）、x（t-2）、……、x（t-k+1）的影响，而这k-1个随机变量又都和x（t-k）具有相关关系，所以自相关系数p（k）里实际掺杂了其他变量对x（t）与x（t-k）的影响。
为了能单纯测度x（t-k）对x（t）的影响，引进偏自相关系数的概念。对于平稳时间序列{x（t）}，所谓滞后k偏自相关系数指在给定中间k-1个随机变量x（t-1）、x（t-2）、……、x（t-k+1）的条件下，或者说，在剔除了中间k-1个随机变量x（t-1）、x（t-2）、……、x（t-k+1）的干扰之后，x（t-k）对x（t）影响的相关程度。
简单来说，就是自相关系数ACF还包含了其他变量的影响，而偏自相关系数PACF是严格这两个变量之间的相关性。在ACF中存在着线性关系和非线性关系，偏自相关函数就是把线性关系从自动关系性中消除。当PACF近似于0，表明两个时间点之间的关系性是完全由线性关系所造成的。
通过R语言的pacf（）函数来进行偏自相关函数计算。
> pacf(tsx)$acf
, , 1
[,1]
[1,] 0.987861891
[2,] 0.006463542
[3,] -0.030541593
[4,] -0.041290415
[5,] 0.047921168
[6,] -0.009774246
[7,] -0.006267004
[8,] 0.002146693
[9,] 0.028782423
[10,] 0.014785187
[11,] 0.019307564
[12,] 0.060879259
[13,] -0.007254278
[14,] -0.004139848
[15,] 0.015707900
[16,] -0.018615370
[17,] 0.037067452
[18,] 0.019322565
[19,] -0.048471479
[20,] 0.023388065
[21,] 0.027640953
[22,] 0.051177900
[23,] 0.028063875
[24,] -0.003957142
[25,] 0.034030631
[26,] 0.004270416
[27,] -0.029613088
[28,] 0.033715973
[29,] 0.092337583
[30,] -0.031264028
# 可视化输出
> pacf(tsx)
从图2-30的结果分析，当滞后为1时AR模型显著，滞后为其他值是PACF的值接近于0不显著。所以，对于数据集tsx来说，数据满足AR（1）的自回归模型。对于上文中参数估计出的1阶自相关系数值是可以用的。
图2-30 PACF检验
通过模型识别，我们已经确定了数据集tsx是符合AR（1）的建模条件的，同时我们也创建了AR（1）模型。接下来，就可以利用这个自回测的模型进行预测，通过规律发现价值。在R语言中，我们可以用predict（）函数实现预测的计算。
使用AR（1）模型进行预测，并保留前5个预测点。
> predict(a,10,n.ahead=5L) $pred Time Series: Start = 2 End = 6 Frequency = 1 [1] 9.839680 9.681307 9.524855 9.370303 9.217627 $se Time Series: Start = 2 End = 6 Frequency = 1 [1] 1.080826 1.519271 1.849506 2.122810 2.359189
上面结果中，变量$pred表示预测值，变量$se为误差。我们可以生成可视化的图，更直观的看到预测的结果。
# 生成50个预测值 > tsp<-predict(a,n.ahead=50L) # 把原数据画图 > plot(tsx) # 把预测值和误差画出来 > lines(tsp$pred,col='red') > lines(tsp$pred+tsp$se,col='blue') > lines(tsp$pred-tsp$se,col='blue')
从图2-31中看到，黑色线为原始数据，在右侧红色线为预测值，蓝色线为预测值的范围。这样我们就利用AR（1）模型，实现了对规律的预测计算。
图2-31 自回归模型预测
上面关于预测和可视化的过程，我们是通过原生的predict（）函数和plot（）函数完成的。在R语言中，可以用forecast包来简化上面的操作过程，让代码更少，操作更便捷。
# 加载forecast包
> library('forecast')
# 生成模型AR(1)
> a2 <- arima(tsx, order=c(1,0,0))
> tsp2<-forecast(a2, h=50)
> plot(tsp2)
图2-32 使用forecast包进行预测可视化
查看forecast（）计算后的预测结果。
> tsp2
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
1001 -15.71590 -16.99118 -14.440628 -17.66627 -13.7655369
1002 -15.60332 -17.39825 -13.808389 -18.34843 -12.8582092
1003 -15.49181 -17.67972 -13.303904 -18.83792 -12.1456966
1004 -15.38136 -17.89579 -12.866932 -19.22685 -11.5358726
1005 -15.27197 -18.06994 -12.474000 -19.55110 -10.9928432
1006 -15.16362 -18.21425 -12.112996 -19.82915 -10.4980922
1007 -15.05631 -18.33593 -11.776682 -20.07206 -10.0405541
1008 -14.95001 -18.43972 -11.460312 -20.28705 -9.6129750
1009 -14.84474 -18.52891 -11.160567 -20.47919 -9.2102846
1010 -14.74046 -18.60591 -10.875013 -20.65216 -8.8287673
1011 -14.63718 -18.67257 -10.601802 -20.80877 -8.4655994
1012 -14.53489 -18.73030 -10.339486 -20.95121 -8.1185723
1013 -14.43357 -18.78024 -10.086905 -21.08123 -7.7859174
1014 -14.33322 -18.82333 -9.843112 -21.20026 -7.4661903
1015 -14.23383 -18.86034 -9.607319 -21.30947 -7.1581923
1016 -14.13538 -18.89190 -9.378864 -21.40985 -6.8609139
通过forecast（）函数，我们得到了直接生成的Forecast值、80%概率的预测值范围和95%概率的预测值范围。
在明白了整个自回归模型的设计思路、建模过程、检验条件、预测计算、可视化展示的完整操作后，我们就可以真正地把自回归模型用到实际的业务中。发现规律，发现价值！
自回归模型只是开始，更多扩展内容请订阅我的博客（http://fens.me）。要把时间序列研究透彻，才能解读出深层的中国经济。