发帖
雷达卡
Lecture3-15
Example Detailed analysis of J&J earnings.
R Demonstration: output edited.
> setwd("D:\\R")
>x=ts(scan("q-earn-jnj.txt"),frequency=4,start=c(1960,1))设置了1960年为基期，后面逐步
Read 84 items
> plot(x)##画图
> y=log(x) ##生成y=lnx
> plot(y) ##画图
> c1=paste(c(1:4))##生成一个新序列，1到4
> points(y,pch=c1) ##给每个点标上1234，顺序循环标法。
> par(mfcol=c(2,1)) ##一页展示两张图。
> acf(y,lag.max=16) ##y的自相关函数，衰减到0非常的慢，4期。
> y1=as.vector(y)##y向量化
> acf(y1,lag.max=16) ##y1的自相关函数，向量化后衰减到0非常的慢，16期
> dy1=diff(y1) ##对向量化的y1进行差分，一阶。常规差分
> acf(dy1,lag.max=16)##差分后看自相关函数。长得像正弦波，
> sdy1=diff(dy1,4) ##对向量化的dy1进行差分，四阶。季节差分
> acf(sdy1,lag.max=12) ##差分后看自相关函数。一阶后截尾。
> m1=arima(y1,order=c(0,1,1),seasonal=list(order=c(0,1,1),period=4))
##建模，常规加季节性。
> m1
Call:
arima(x = y1, order = c(0, 1, 1), seasonal= list(order = c(0, 1, 1), period = 4))
Coefficients:
ma1 sma1
-0.6809 -0.3146
s.e. 0.0982 0.1070
PS:
(1-B^4)(1-B)R(t)=(1-0.3146^4)(1-0.6809B)a(t)
Var[a(t)]=0.00793
sigma^2 estimated as 0.007931: log likelihood = 78.38, aic = -150.75
> par(mfcol=c(1,1))##一页一图
> tsdiag(m1)##模型检测图
> f1=predict(m1,8)##预测
> names(f1)
[1] "pred" "se"
> f1
$pred
Time Series:
Start = 85
End = 92
Frequency = 1
[1] 2.905343 2.823891 2.912148 2.5810853.036450 2.954999 3.043255 2.712193
$se
Time Series:
Start = 85
End = 92
Frequency = 1
[1] 0.08905414 0.09347899 0.097703660.10175307 0.13548771 0.14370561
[7] 0.15147833 0.15887123
Example.U.S.weekly interest rate data:
1-year and 3-year Constant maturity rates.
Data are shown in Figure 6.
R Demonstration: output edited.
> setwd("D:\\R")##设置目录
>da=read.table("w-gs1n36299.txt",header=T)##读取数据，记住读行标题，因为原数据带标题，不读取标题，后续都会报错，相当于数据被引入了部分文本的内容。
> r1=da[,1]##生成 数据
> r3=da[,2]##生成数据
> lines(1:1967,r3,lty=2)##没跑出来，不明白为什么
Error in plot.xy(xy.coords(x, y), type =type, ...) :
plot.new has not been called yet
> m1=lm(r3~r1)##直接跑回归模型，
> summary(m1)
Call:
lm(formula = r3 ~ r1)
Residuals:
Min 1Q Median 3Q Max
-1.8121 -0.4023 0.0031 0.4026 1.3388
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.910687 0.032250 28.24 <2e-16 ***
r1 0.923854 0.004389 210.51 <2e-16 ***
r3=0.91+0.92r1
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’1
Residual standard error: 0.538 on 1965degrees of freedom
Multiple R-squared: 0.9575, Adjusted R-squared: 0.9575
F-statistic: 4.431e+04 on 1 and 1965DF, p-value: < 2.2e-16
> acf(m1$residual)##残差自相关图，30期之后仍然没有衰减到0，很慢。
> c3=diff(r3)##
> c1=diff(r1)##分别对该序列进行差分
> par(mfcol=c(2,1))
> plot(r1,r3)
> plot(c1,c3)##输出在一页上方便比较。
> m2=lm(c3~c1)##用c3和c1直接建模
> summary(m2)###输出汇总结果
Call:
lm(formula = c3 ~ c1)
Residuals:
Min 1Q Median 3Q Max
-0.38060 -0.03338 -0.00054 0.03437 0.47418
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.0002475 0.0015380 0.161 0.872
c1 0.7810590 0.0074651 104.628 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’1
PS: c3=0.0002475+0.781c1
Residual standard error: 0.06819 on 1964degrees of freedom
Multiple R-squared: 0.8479, Adjusted R-squared: 0.8478
F-statistic: 1.095e+04 on 1 and 1964DF, p-value: < 2.2e-16
> acf(m2$residual)
> acf(m1$residual)##比较两者的残差的自相关图，改善明显，一阶截尾。
> plot(m2$residuals,type="l")
> m3=arima(c3,xreg=c1,order=c(0,0,1))，残差用MA1来拟合。
> m3
Call:
arima(x = c3, order = c(0, 0, 1), xreg =c1)
Coefficients:
ma1 intercept c1
0.2115 0.0002 0.7824
s.e. 0.0224 0.0018 0.0077
PS: c3=0.0002+0.7824c1+a(t)+0.2115a(t-1)
sigma^2 estimated as 0.004456: log likelihood = 2531.84, aic = -5055.69
> acf(m3$residual)##残差相关图
> tsdiag(m3)##模型检验图
> m4=arima(c3,xreg=c1,order=c(1,0,0))##残差用AR1来估计
> m4
Call:
arima(x = c3, order = c(1, 0, 0), xreg =c1)
Coefficients:
ar1 intercept c1
0.1922 0.0003 0.7829
s.e. 0.0221 0.0019 0.0077
PS: c3=0.0003+0.7829c1+a(t)+0.1922a(t-1)+e(t)
其中：a(t)=0.1922a(t-1)+e(t)
从AIC看还是MA模型更优一些。
sigma^2 estimated as 0.004474: log likelihood = 2527.86, aic = -5047.72
> acf(m3$residual)
> acf(m4$residual)##比较自相关图
> tsdiag(m4) ##模型检验图
提升卡
置顶卡
沉默卡
变色卡
抢沙发
千斤顶
显身卡
京公网安备 11010802022788号
