时间序列分析——基于R（第2版）—第7章

7-1 某地区过去38年的谷物产量（y）和这些年该地区相应的降雨量（x）

1. 使用单位根检验分别考察这两个模型的平稳性
2. 选择适当模型分别拟合这两个序列的发展
3. 确定这两个序列之间是否具有协整关系
4. 如果这两个序列之间具有协整关系，建立适当的模型拟合谷物产量序列发展

x=ts(E7_1$x)
y=ts(E7_1$y)
plot(x)
plot(y)
library(aTSA)
adf.test(x)#adf检验的原假设是存在单位根，序列不平稳
#p<0.05,拒绝原假设，序列平稳
for(i in 1:6) print(Box.test(x,lag=3*i))
#纯随机性检验 p值>0.05,接受原假设，属于白噪声序列
adf.test(y)#平稳
for(i in 1:6) print(Box.test(y,lag=3*i))
#两个序列皆平稳
#两个序列皆为白噪声序列
#两变量都具有白噪声属性，所以不考虑协整方程

7-2 在一定浓度溶液中，考察草履虫和某种草履虫掠食动物之间动态数量的变化

1. 考察这两种生物之间的动态关系，检验他们之间是否具有协整关系
2. 选择适当的模型拟合这两种生物之间的动态关系，并预测未来一周这两种生物的浓度

library(forecast)
library(aTSA)
x=ts(E7_2$x)
y=ts(E7_2$y)
plot(x)
lines(y,col=2)
coint.test(y,x)#协整检验,P<0.01,具有协整关系
x1=x[-c(1,2)]
y1=y[-c(37,38)]
fit=Arima(y1,xreg = x1,order=c(2,0,0))
fit
tsdiag(fit)
acf(diff(x,10))
pacf(diff(x,10))
fitx=Arima(x,order = c(2,0,0),seasonal = list(order=c(0,1,1),period=10))
fitx
tsdiag(fitx)
forex=forecast::forecast(fitx,h=14);forex
plot(forex)
lines(forex$fitted,col=2)
forey=forecast::forecast(fit,xreg=forex$mean,h=14);forey
plot(forey)
lines(forey$fitted,col=2)
#具有协整关系
#掠食动物和被掠食动物都显示出显著的周期特征两个序列皆为非平稳序列；但是掠食动物和被掠食动物延迟2阶具有协整关系

7-3我国1950-2008年进出口总额

1. 使用单位根检验分别考察这进口总额和出口总额序列的平稳性
2. 分别对进口总额序列和出口总额序列拟合模型
3. 考察这两个序列之间是否具有协整关系
4. 如果这两个序列之间具有协整关系，建立适当的模型拟合相关关系
5. 构造该协整模型的误差修正模型

library(forecast)
library(aTSA)
x=ts(E7_3$x,start = 1950)
y=ts(E7_3$y,start = 1950)
plot(x)
lines(y,col=2)
acf(diff(diff(x)))
pacf(diff(diff(x)))
fitx=arima(x,order=c(0,2,0))#进口
fitx
ts.diag(fitx)
acf(diff(diff(y)))
pacf(diff(diff(y)))
fity=arima(y,order=c(2,2,0),transform.pars = F,fixed = c(0,NA))#出口
fity
ts.diag(fity)
coint.test(y,x)#协整检验
fit=Arima(y,xreg = x,order=c(2,0,0))
fit
tsdiag(fit)
ecm(y,x)#构造该协整模型的误差修正模型

7-4我国1979-2014年社会消费品零售总额序列和国内生产总值序列数据

1. 分别对这两个序列拟合arima模型，并预测未来5年的序列发展
2.  考察这两个序列之间是否具有协整关系
3.  如果存在协整关系，请思考这两个序列之间的因果关系（哪个是自变量，哪个是因变量）,并构造协整模型，预测未来5年的序列发展

library(forecast)
library(aTSA)
x=ts(E7_4$x,start = 1979)
y=ts(E7_4$y,start = 1979)
plot(x)
plot(y)
acf(diff(diff(x)))
pacf(diff(diff(x)))
fitx=arima(x,order=c(3,2,0),transform.pars = F,fixed = c(0,0,NA))
fitx
ts.diag(fitx)
forex=forecast::forecast(fitx,h=5);forex
acf(diff(diff(y)))
pacf(diff(diff(y)))
fity=arima(y,order=c((2),2,0),transform.pars = F,fixed = c(0,NA));fity
ts.diag(fity)
forey=forecast::forecast(fity,h=5);forey
coint.test(y,x)
#不存在协整关系

7-5 英国1969 1月-1984年12月车祸数据

1. 研究安全带强制法律的执行是否对司机伤亡数据有显著的干预作用。
2. 研究司机伤亡数据与行驶里程、汽油价格及安全带强制法律的执行之间是否具有协整关系．
3. 研究安全带强制法律的执行是否对前座乘客伤亡数据有显著的干预作用。
4. 研究前座乘客伤亡数据与行驶里程、汽油价格及安全带强制法律的执行之间是否具有协整关系．
5. 研究安全带强制法律的执行是否对后座乘客伤亡数据有显著的干预作用．meiyou
6. 研究后座乘客伤亡数据与行驶里程、汽油价格及安全带强制法律的执行之间是否具有协整关系．
7. 研究司机伤亡数据、前座乘客伤亡数据和后座乘客伤亡数据之间是否具有协整关系

library(forecast)
library(aTSA)
Seatbelts#z自带数据集
drivers=Seatbelts[,2]
front=Seatbelts[,3]
rear=Seatbelts[,4]
kms=Seatbelts[,5]
petroprice=Seatbelts[,6]
law=Seatbelts[,8]
x=matrix(Seatbelts[,c(5,6,8)],ncol=3)
coint.test(drivers,x)
fit1=arima(drivers,xreg = data.frame(petroprice,law),order=c(2,0,0),
           seasonal=list(order=c(1,1,0),period=12))
fit1
ts.diag(fit1)
#安全带强制法律的执行对司机伤亡数据有显著的干预作用
#司机伤亡数据与行驶里程、汽油价格及安全带强制法律的执行之间具有协整关系
coint.test(front,x)
fit2=arima(front,xreg = data.frame(petroprice,law),order=c(2,0,0),
           seasonal = list(order=c(1,1,0),period=12))           
fit2
ts.diag(fit2)
#安全带强制法律的执行对前座乘客伤亡数据有显著的干预作用
#前座乘客伤亡数据与行驶里程、汽油价格及安全带强制法律的执行之间具有协整关系．
coint.test(rear,x)
fit3=arima(rear,xreg = data.frame(kms,petroprice),order=c(1,0,0),
           seasonal = list(order=c(1,1,0),period=12))
fit3
ts.diag(fit3)

y=matrix(Seatbelts[,c(3,4)],ncol=2)
coint.test(drivers,y)

7-6 1867-1947年玉米价格、玉米产量、生猪价格、生猪产量以及农场工人平均工资的数据，

1. 分析这五个变量的单整情况．
2. 分析这五个变量的 Granger 因果关系．
3. 分析农场工人平均工资、农场作物的价格、家畜的价格与供应量之间是否具有协整关系．如果有，拟合协整模型与误差修正模型，并解释这两个模型中各参数的意义.

library(forecast)
library(aTSA)
x1=ts(E7_6$maize_price,start=1967)
x2=ts(E7_6$maize_yield,start=1967)
x3=ts(E7_6$pig_price,start=1967)
x4=ts(E7_6$pig_yield,start=1967)
x5=ts(E7_6$salary,start=1967)
plot(x1)
adf.test(x1)
adf.test(diff(x1))
plot(x2)
adf.test(x2)
adf.test(diff(x2))
plot(x3)
adf.test(x3)
adf.test(diff(x3))
plot(x4)
adf.test(x4)
adf.test(diff(x4))
plot(x5)
adf.test(x5)
adf.test(diff(x5))
x=data.frame(x1,x2,x3,x4)
x=as.matrix(x)
coint.test(x5,x,d=1)
fit=arima(x5,xreg = data.frame(x1,x3,x4),order = c(3,0,0),transform.pars = F,
          fixed = c(NA,0,NA,NA,NA,NA,NA))
fit
ts.diag(fit)
ecm(x5,as.matrix((data.frame(x1,x3,x4))))
#(1）这五个变量均为1阶单整序列
#(2）农场工人薪水为果，玉米价格、玉米产量、生猪价格、生猪产量为因。
#(3）分析农场工人平均工资、农场作物、家畜的价格与供应链之间具有协整关系。拟合协整模型如fit

7-7

1. 分别绘制这四个序列的时序图，考察这四个序列各自的波动特征，研究它们的单整性，并分别拟合单变量 ARIMA 模型．
2. 考察这四个变量的 Granger 因果关系．
3. 以 GNP 为响应序列，根据因果检验结果选择适当的自变量，考察自变量与响应变量之间是否具有协整关系
4. 如果这些宏观经济变量之间具有协整关系，则拟合协整模型与误差修正模型并解释这两个模型的系数意义

library(forecast)
library(aTSA)
logM1=ts(E7_7$logM1,start = c(1954,4),frequency = 4)
logGNP=ts(E7_7$logGNP,start = c(1954,4),frequency = 4)
sr=ts(E7_7$sr,start = c(1954,4),frequency = 4)
lr=ts(E7_7$lr,start = c(1954,4),frequency = 4)
plot(logM1)
adf.test(logM1)
adf.test(diff(logM1))
acf(diff(logM1))
pacf(diff(logM1))
fit1=arima(logM1,order = c(1,1,0))
fit1
ts.diag(fit1)

plot(logGNP)
adf.test(logGNP)
adf.test(diff(logGNP))
acf(diff(logGNP))
pacf(diff(logGNP))
fit2=arima(logGNP,order = c(2,1,1))
fit2
ts.diag(fit2)

plot(sr)
adf.test(sr)
adf.test(diff(sr))
acf(diff(sr))
pacf(diff(sr))
fit3=arima(sr,order = c(7,1,1),transform.pars = F,fixed = c(0,NA,0,0,0,0,NA,NA))
fit3
ts.diag(fit3)

plot(lr)
adf.test(lr)
adf.test(diff(lr))
acf(diff(lr))
pacf(diff(lr))
fit4=arima(lr,order = c(14,1,0),transform.pars = F,fixed = c(NA,0,0,0,0,0,0,0,0,0,0,0,0,NA))
fit4
ts.diag(fit4)

lagGNP=lag(logGNP)
y=logGNP[1:134]
x=data.frame(logM1[1:134],sr[1:134],lr[1:134],lagGNP[2:135])
x=as.matrix(x)
coint.test(y,x)
fit=arima(y,xreg = data.frame(sr[1:134],lr[1:134],lagGNP[2:135]))
fit
ts.diag(fit)
ecm(y,as.matrix(data.frame(sr[1:134],lr[1:134],lagGNP[2:135])))
#分别绘制这四个序列的时序图，这四个变量均为1阶单整序列
#考察这四个变量的 Granger 因果关系．logGNP为果，其他三个序列为因
#协整模型为 0.0683+0.9893logGNPt-1+0.3257srt-0.1650lrt
#误差修正模型为：0.50876sr+0.13727lr+0.97062logGNPt-1-0.89442ECMt-1