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

5-1 我国1949-2008年每年铁路货运量

1. 选择适当的模型拟合该序列，并预测2009-2013年我国铁路货运量

x=ts(E5_1$volume,start = 1949)
plot(x)
plot(x)
plot(diff(x))
library(tseries)
adf.test(diff(x))#tseries包
Box.test(diff(x))
acf(diff(x))
pacf(diff(x))

#该序列1阶差分后平稳，根据差分后序列的自相关图拖尾和偏自相关图1阶截尾特征，对该序列拟合 ARIMA (1,1,0）模型
fit=Arima(x,order=c(1,1,0),include.drift=T)
fit
tsdiag(fit)
library(forecast)
fore=forecast(fit,h=5);fore
par(mfrow=c(1,1))
plot(fore)
lines(fore$fitted,col=2)

5-2 1750-1849年瑞典人口出生率数据
 

x=ts(E5_2$birth_rate,start=1750)
plot(x)
adf.test(x)
Box.test(x)
library(tseries)
acf(x)
pacf(x)

#该序列 adf 检验平稳，根据自相关图和偏自相关图特征，可以识别为偏自相关图1阶截尾，拟合 AR (1）模型
fit=arima(x,order=c(1,0,0));fit

5-3 1867-1938年英国的绵羊数量

1. 确定该序列的平稳性
2. 选择适当的模型拟合该序列的发展
3. 利用拟合模型预测1939-1945年英国的绵羊数量

x=ts(E5_3$number,start=1867)
plot(x)
plot(diff(x))
library(tseries)
adf.test(diff(x))#tseries
Box.test(diff(x))
acf(diff(x))
pacf(diff(x))

#1根据时序图和 adf 检验结果可知，原序列非平稳，一阶差分后序列平稳
library(forecast)
fit1=Arima(x,order = c(3,1,0),include.drift=T);fit1

#2一阶差分后序列自相关图拖尾，偏自相关图3阶截尾，所以对原序列拟合模型 ARIMA (3,1,0)，参数输出结果显示2阶自相关系数与均值都不显著非零。所以最后拟合不带均值项的疏系数模型 ARIMA ((1,3),1,0)

fit2=arima(x,order = c(3,1,0),transform.pars = F,fixed=c(NA,0,NA));fit2
tsdiag(fit2)

#3 1939-1945年英国的绵羊数量预测结果如下
fore=forecast(fit2,h=7);fore
plot(fore)
lines(fore$fitted,col=2)

5-4我国人口出生率，死亡率和自然增长率

1. 分析我国人口出生率，死亡率和自然增长率的平稳性
2. 对非平稳序列选择适当的差分方式实现差分后平稳
3. 选择适当的模型拟合我人口出生率的变化，并预测未来10年的人口出生率。
4. 选择适当的模型拟合我国人口死亡率的变化，并预测未来10年的人口死亡率
5. 选择适当的模型拟合我国人口自然增长率的变化，并预测未来10年的人口自然增长率

  

x1=ts(E5_4$birth_rate,start = 1980)
x2=ts(E5_4$mortality,start = 1980)
x3=ts(E5_4$ngr,start = 1980)
plot(x1)
plot(diff(x1))
library(tseries)
adf.test(diff(x1))
for(i in 1:6) print(Box.test(diff(x1),i))
acf(diff(x1))
pacf(diff(x1))
fit11=arima(x1,order=c(0,1,0))
fit12=arima(x1,order=c(0,1,1))
fit13=arima(x1,order=c(1,1,0))
c(AIC(fit11),AIC(fit12),AIC(fit13),BIC(fit11),BIC(fit12),BIC(fit13))
fit12
library(forecast)
fore1=forecast(fit12,h=10);fore1

plot(x2)
plot(diff(x2))
library(tseries)
adf.test(diff(x2))
for(i in 1:6) print(Box.test(diff(x2),i))
acf(diff(x2))
pacf(diff(x2))
fit21=arima(x2,order=c(0,1,0))
fit22=arima(x2,order=c(0,1,1))
fit23=arima(x2,order=c(1,1,0))
c(AIC(fit21),AIC(fit22),AIC(fit23),BIC(fit21),BIC(fit22),BIC(fit23))
fit21
library(forecast)
fore2=forecast(fit21,h=10);fore2

plot(x3)
plot(diff(x3))
library(tseries)
adf.test(diff(x3))
for(i in 1:6) print(Box.test(diff(x3),i))
acf(diff(x3))
pacf(diff(x3))
fit31=arima(x3,order=c(0,1,0))
fit32=arima(x3,order=c(0,1,1))
fit33=arima(x3,order=c(1,1,0))
c(AIC(fit31),AIC(fit32),AIC(fit33),BIC(fit31),BIC(fit32),BIC(fit33))
fit32
library(forecast)
fore3=forecast(fit32,h=10);fore3

#1我国人口出生率、死亡率和自然增长率序列均有显著线性趋势，为典型非平稳序列。
#2根据图检验和 adf 检验结果显示，我国人口出生率、死亡率和自然增长率序列的一阶差分后序列均平稳。
#3我国人口出生率可以拟合 ARIMA (0,1,0）模型，也可以拟合 ARIMA (0,1,1）模型。其中 ARIMA (0,1,1）模型 AIC , BIC 更小一点。得到的拟合模型参数为。基于该模型预测值为
#4我国人口死亡率可以拟合 ARIMA (0,1,0）模型，预测，，
#5我国人口自然增长率可以拟合 ARIMA (0,1,0）模型， ARIMA (1,1,0）模型，也可以拟合 ARIMA (0,1,1）模型。其中 ARIMA (0,1,1）模型 AIC , BIC更小一点。得到的拟合模型参数为，基于 ARIMA (0.1,1）模型，预测1998年之后人口自然增长率未来十年的预测值为
                                      

5-5 某农场1867-1947年玉米和生猪的销售价格、产量及农场工人平均收人如表

1. 分析这几个序列的平稳性．
2. 对非平稳序列找到适当的差分阶数实现差分后平稳．
3. 选择适当的模型拟合这几个序列的发展，并做10年期序列预测．
4. 绘制拟合与预测效果图

x1=ts(E5_5$maize_price,start = 1867)
x2=ts(E5_5$maize_yield,start = 1867)
x3=ts(E5_5$wages,start = 1867)
x4=ts(E5_5$pig_price,start = 1867)
x5=ts(E5_5$pig_yield,start = 1867)

plot(x1)
plot(diff(x1))
for(i in 1:6) print(Box.test(diff(x1),i))
acf(diff(x1))
pacf(diff(x1))
#对玉米价格序列拟合疏系数模型 ARIMA (0.1,(7))，参数估计值如下
fit1=arima(x1,order=c(0,1,7),transform.pars = F,fixed = c(0,0,0,0,0,0,NA))
fit1
library(forecast)
fore1=forecast(fit1,h=10);fore1#对玉米价格序列进行10期预测，预测值如下：
plot((fore1))#预测效果图

plot(x2)
plot(diff(x2))
for(i in 1:6) print(Box.test(diff(x2),i))
acf(diff(x2))
pacf(diff(x2))
#对玉米产量拟合arima（0，1，1）模型
fit2=arima(x2,order=c(0,1,1));fit2
library(forecast)
fore2=forecast(fit2,h=10);fore2
plot((fore2))#预测效果图

plot(x3)
plot(diff(x3))
for(i in 1:6) print(Box.test(diff(x3),i))
acf(diff(x3))
pacf(diff(x3))
#对工人工资序列拟合疏系数模型 ARIMA ((1,13),1,0)，参数估计值如下
fit3=arima(x3,order=c(13,1,0),transform.pars = F,fixed = c(NA,0,0,0,0,0,0,0,0,0,0,0,NA))
fit3
library(forecast)
fore3=forecast(fit3,h=10);fore3
plot((fore3))#预测效果图

plot(x4)
plot(diff(x4))
for(i in 1:6) print(Box.test(diff(x4),i))
acf(diff(x4))
pacf(diff(x4))
#对生猪价格拟合arima（3，1，0）模型
fit4=arima(x4,order=c(3,1,0));fit4
library(forecast)
fore4=forecast(fit4,h=10);fore4
plot((fore4))#预测效果图

plot(x5)
plot(diff(x5))
for(i in 1:6) print(Box.test(diff(x5),i))
acf(diff(x5))
pacf(diff(x5))
#对生猪产量拟合疏系数模型arima（（2，3，5），1，0）
fit5=arima(x5,order=c(5,1,0),transform.pars = F,fixed = c(0,NA,NA,0,NA))
fit5
library(forecast)
fore5=forecast(fit5,h=10);fore5
plot((fore5))#预测效果图

#1玉米价格、玉米产量、工人薪水、生猪价格、生猪产量序列都不平稳，
#2玉米价格、玉米产量、工人薪水、生猪价格、生猪产量序列一阶差分都实现平稳。