金融时间序列得线性模型——自回归.doc

金融时间序列的线性模型——自回归R实例例2.3> setwd("C:/Users/Mr.Cheng/Desktop/课件/金融数据分析导论基于R/DataSets/ch2data")%设置工作> da=read.table("q-gnp4710.txt",header=T)> head(da) Year Mon Dat VALUE1 1947 1 1 238.12 1947 4 1 241.53 1947 7 1 245.64 1947 10 1 255.65 1948 1 1 261.76 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)> m1Call:arima(x = gnp, order = c(3, 0, 0))Coefficients: ar1 ar2 ar3 intercept 0.4386 0.2063 -0.1559 0.0163s.e. 0.0620 0.0666 0.0626 0.0012sigma^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.864212i> 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" "aic" [6] "n.used" "order.max" "partialacf" "resid" "method" > 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.09290084l 模型的检验> vw=read.table('m-ibm3dx2608.txt',header=T)[,3]> m3=arima(vw,order=c(3,0,0))%用AR〔3〕拟合> m3Call:arima(x = vw, order = c(3, 0, 0))Coefficients: ar1 ar2 ar3 intercept 0.1158 -0.0187 -0.1042 0.0089s.e. 0.0315 0.0317 0.0317 0.0017sigma^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 testdata: m3$residualsX-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> m3Call: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.0089s.e. 0.0313 0 0.0315 0.0017sigma^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 testdata: m3$residualsX-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 length> tratio ar1 ar2 ar3 intercept 3.6301072 0.0000000 -62.0713895 0.2859641 画自相关函数> 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)。