金融时间序列的线性模型——自回归R实例
例2.3
> setwd("C:/Users/Mr.Cheng/Desktop/课件/金融数据分析导论基于R/Dat
aSets/ch2data")%设置工作目录
> da=read.table("q-gnp4710.txt",header=T)
> head(da)
Year Mon Dat VALUE
1 1947 1 1 238.1
2 1947 4 1 241.5
3 1947 7 1 245.6
4 1947 10 1 255.6
5 1948 1 1 261.7
6 1948 4 1 268.7
> G=da$VALUE
> LG=log(G)
> gnp=diff(LG)
> dim(da)
[1] 253 4
> tdx=c(1:253)/4+1947 %创建一个时间序列指数，从1947开始，每次增加一个季度，一共253个季度。
> par(mfcol=c(2,1))画两行一列的小图
> plot(tdx,LG,xlab='year',ylab='GNP',type="l
> plot(tdx[2:253],gnp,type='l',xlab='year',ylab='growth')
> acf(gnp,lag=12)%画滞后12阶的对数增长率的自相关图
> pacf(gnp,lag=12)%画滞后12阶的对数增长率的偏自相关图
> m1=arima(gnp,order=c(3,0,0))%计算AR(3)
> m1
Call:
arima(x = gnp, order = c(3, 0, 0))
Coefficients:
ar1 ar2 ar3 intercept
0.4386 0.2063 -0.1559 0.0163
s.e. 0.0620 0.0666 0.0626 0.0012
sigma^2 estimated as 9.549e-05: log likelihood = 808.56, aic = -1607.12
> tsdiag(m1,gof=12)%模型检验
> p1=c(1,-m1$coef[1:3])%设置多项式方程的系数：
1-0.438z-0.206z2+0.156z3=0
> r1=polyroot(p1)%解多项式方程得到特征根
> r1
[1] 1.616116+0.864212i -1.909216-0.000000i 1.616116-0.864212 i
> Mod(r1)%计算特征根的模
[1] 1.832674 1.909216 1.832674
> k=2*pi/acos(1.616116/1.832674)%计算周期
> k
[1] 12.79523
> mm1=ar(gnp,method='mle')%用AIC准则自动为AR（P）定阶，方法为极大似然估计
> mm1$order%查看阶数
[1] 9
> names(mm1)%得到mm1的名字
[1] "order" "ar" "var.pred" "x.mean" "a ic"
[6] "https://www.doczj.com/doc/7d10314345.html,ed" "order.max" "partialacf" "resid" " method"
[11] "series" "frequency" "call" "asy.var.coef" > print(mm1$aic,digits = 3)%查看mm1中的aic值，保留三位小数
0 1 2 3 4 5 6 7 8 9 10 11
77.767 11.915 8.792 4.669 6.265 5.950 5.101 4.596 6.541 0.000 0.509 2.504
12
2.057
> aic=mm1$aic
> length(aic)
[1] 13
> plot(c(0:12),aic,type='h',xlab='order',ylab='aic')%画aic竖线
图
> lines(0:12,aic,lty=2)%画aic连线图（虚线）
> vw=read.table('m-ibm3dx2608.txt',header=T)[,3]%读取第3列数据
> t1=prod(vw+1)%计算35年后的终值
> t1
[1] 1592.953
> head(vw)
[1] 0.000724 -0.033374 -0.064341 0.038358 0.012172 0.056888
> t1^(12/996)-1%折算回平均每年的回报
[1] 0.09290084
模型的检验
> vw=read.table('m-ibm3dx2608.txt',header=T)[,3]
> m3=arima(vw,order=c(3,0,0))%用AR（3）拟合
> m3
Call:
arima(x = vw, order = c(3, 0, 0))
Coefficients:
ar1 ar2 ar3 intercept
0.1158 -0.0187 -0.1042 0.0089
s.e. 0.0315 0.0317 0.0317 0.0017
sigma^2 estimated as 0.002875: log likelihood = 1500.86, aic = -2991.73
> (1-.1158+.0187+.1042)*mean(vw)%计算phi(0)
[1] 0.008967611
> sqrt(m3$sigma2)%计算残差标准误
[1] 0.0536189
> Box.test(m3$residuals,lag=12,type="Ljung")%检验残差的自相关函数，如果显示出额外的序列相关性，则应该考虑到这些相关性并进行扩展Box-Ljung test
data: m3$residuals
X-squared = 16.352, df = 12, p-value = 0.1756
> pv=1-pchisq(16.35,9)%由上一步算得Q（12）=16.352,并且基于它所渐进服从的自由度为9(修正自由度12-2)的卡方分布，得到p值为0.06，因此在5%的显著水平下无法拒绝原假设
> pv
[1] 0.05992276
> m3=arima(vw,order=c(3,0,0),fixed=c(NA,0,NA,NA))%改进模型：由于间隔为2的AR系数在5%的水平下不显著，因此修改后的模型去除2阶滞后项。（下面有补充计算）
Warning message:
In arima(vw, order = c(3, 0, 0), fixed = c(NA, 0, NA, NA)) :
一些AR参数是固定的：把transform.pars设成FALSE
> m3
Call:
arima(x = vw, order = c(3, 0, 0), fixed = c(NA, 0, NA, NA)) Coefficients:
ar1 ar2 ar3 intercept
0.1136 0 -0.1063 0.0089
s.e. 0.0313 0 0.0315 0.0017
sigma^2 estimated as 0.002876: log likelihood = 1500.69, aic = -2993.38
> (1-.1136+.1063)*.0089 %计算phi(0)
[1] 0.00883503
> sqrt(m3$sigma2)
[1] 0.05362832
> Box.test(m3$residuals,lag=12,type='Ljung')
Box-Ljung test
data: m3$residuals
X-squared = 16.828, df = 12, p-value = 0.1562
> pv=1-pchisq(16.83,10)%修正自由度（12-2）
> pv
[1] 0.07821131
%改进后的模型对数据的动态线性相依性的建模是充分的。
关于系数显著性的计算：
> vw=read.table('m-ibm3dx2608.txt',header=T)[,3]
> m3=arima(vw,order=c(3,0,0),fixed=c(NA,0,NA,NA))
Warning message:
In arima(vw, order = c(3, 0, 0), fixed = c(NA, 0, NA, NA)) :
一些AR参数是固定的：把transform.pars设成FALSE
> names(m3)
[1] "coef" "sigma2" "var.coef" "mask" "loglik" " aic"
[7] "arma" "residuals" "call" "series" "code" " n.cond"
[13] "nobs" "model"
> tratio=m3$coef/sqrt(diag(m3$var.coef))%diag函数用于提取对角线上
的元素。
Warning message:
In m3$coef/sqrt(diag(m3$var.coef)) :
longer object length is not a multiple of shorter object lengt h
> tratio
ar1 ar2 ar3 intercept
3.6301072 0.0000000 -62.0713895 0.2859641
显著性取0.05时就把|t|和1.96（查正态分布表的0.975对应的值）比较，大于
就显著，小于就不显著。显著性取0.01时对比2.575，显著性取0.1时对比1.64 5.
画自相关函数
> po=1
> p1=0.8
> T=5000
> x=rep(0,T)%重复产生T个0的向量存储在x中。
> a=rnorm(T)
> for(i in 2:T)
+ x[i]=po+p1*x[i-1]+a[i]
> p2=-.8
> y=rep(0:T)
> for(i in 2:T)
+ y[i]=po+p2*y[i-1]+a[i]
> par(mfcol=c(1,2))
> acf(x,lag=12)
> acf(y,lag=12)