1 目的
为研究国民生产总值与货币供应量及利率的关系。现收集到1954年1月至1987年10月M1货币量对数序列log(M1)，美国月度国民生产总值对数序列l log(GNP)，以及短期利率和长期利率序列，部分数据见表1。该篇文章主要对四个序列分别研究其单整性及按照实际情况拟合A R I M A ARIMAARIMA模型，后续将研究研究这些变量的格兰杰因果关系及协整关系。
表1 部分数据展示
2 原序列可视化
运行程序：
rm(list=ls()) #清空变量 library("openxlsx") #加载包 library("knitr") #加载包 library("xlsx") #加载包 library("tseries") #加载包 data<-read.xlsx("F:\\时间序列分析\\数据.xlsx",'Sheet1', encoding ="UTF-8")#读取数据 M1<-ts(data$log.M1.,start = c(1954,1),frequency = 4) #log（M1） GNP<-ts(data$log.GNP.,start = c(1954,1),frequency = 4) #log（GNP） SF<-ts(data$短期利率,start = c(1954,1),frequency = 4) #短期利率 LF<-ts(data$长期利率,start =c(1954,1),frequency = 4) #长期利率 par(mfrow=c(2,2)) #创建2×2画布 plot(M1,sub='图1 log（M1）时序图',ylab="log（M1）") plot(GNP,sub='图2 log（GNP）时序图',ylab="log（GNP）") plot(SF,sub='图3 短期利率时序图') plot(LF,sub='图4 长期利率时序图')
运行结果：
通过各变量的时序图(图1-图4）可以发现，各时序数据具有明显趋势效应。
3 各序列单整性检验
3.1 log（M1）单整性检验
1.平稳性检验
运行程序：
library(aTSA) adf.test(M1) #log（M1）序列平稳性检验
运行结果：
## Augmented Dickey-Fuller Test ## alternative: stationary ## ## Type 1: no drift no trend ## lag ADF p.value ## [1,] 0 2.290 0.990 ## [2,] 1 1.062 0.921 ## [3,] 2 0.899 0.900 ## [4,] 3 0.738 0.854 ## [5,] 4 0.742 0.855 ## Type 2: with drift no trend ## lag ADF p.value ## [1,] 0 1.084 0.990 ## [2,] 1 -0.658 0.816 ## [3,] 2 -0.971 0.707 ## [4,] 3 -1.395 0.558 ## [5,] 4 -1.425 0.547 ## Type 3: with drift and trend ## lag ADF p.value ## [1,] 0 0.324 0.990 ## [2,] 1 -1.291 0.871 ## [3,] 2 -1.610 0.737 ## [4,] 3 -2.057 0.548 ## [5,] 4 -2.082 0.538 ## ---- ## Note: in fact, p.value = 0.01 means p.value <= 0.01
2.一阶差分后平稳性检验
运行程序：
adf.test(diff(M1)) #一阶差分log（M1）序列平稳性检验
运行结果：
## Augmented Dickey-Fuller Test ## alternative: stationary ## ## Type 1: no drift no trend ## lag ADF p.value ## [1,] 0 -5.96 0.01 ## [2,] 1 -4.76 0.01 ## [3,] 2 -3.88 0.01 ## [4,] 3 -3.65 0.01 ## [5,] 4 -3.44 0.01 ## Type 2: with drift no trend ## lag ADF p.value ## [1,] 0 -6.05 0.01 ## [2,] 1 -4.84 0.01 ## [3,] 2 -3.94 0.01 ## [4,] 3 -3.72 0.01 ## [5,] 4 -3.52 0.01 ## Type 3: with drift and trend ## lag ADF p.value ## [1,] 0 -6.14 0.0100 ## [2,] 1 -4.94 0.0100 ## [3,] 2 -4.05 0.0100 ## [4,] 3 -3.84 0.0192 ## [5,] 4 -3.63 0.0328 ## ---- ## Note: in fact, p.value = 0.01 means p.value <= 0.01
3. log（GNP）单整性检验
1.平稳性检验
运行程序：
adf.test(GNP) #log（GNP）序列平稳性检验
运行结果：
## Augmented Dickey-Fuller Test ## alternative: stationary ## ## Type 1: no drift no trend ## lag ADF p.value ## [1,] 0 8.62 0.99 ## [2,] 1 5.13 0.99 ## [3,] 2 4.11 0.99 ## [4,] 3 4.21 0.99 ## [5,] 4 4.05 0.99 ## Type 2: with drift no trend ## lag ADF p.value ## [1,] 0 -0.695 0.803 ## [2,] 1 -0.793 0.769 ## [3,] 2 -0.656 0.817 ## [4,] 3 -0.564 0.849 ## [5,] 4 -0.423 0.899 ## Type 3: with drift and trend ## lag ADF p.value ## [1,] 0 -1.53 0.772 ## [2,] 1 -2.16 0.503 ## [3,] 2 -2.40 0.405 ## [4,] 3 -2.15 0.511 ## [5,] 4 -2.05 0.552 ## ---- ## Note: in fact, p.value = 0.01 means p.value <= 0.01
2.一阶差分后平稳性检验
运行程序：
adf.test(diff(GNP)) #一阶差分log（GNP）序列平稳性检验
运行结果：
## Augmented Dickey-Fuller Test ## alternative: stationary ## ## Type 1: no drift no trend ## lag ADF p.value ## [1,] 0 -6.13 0.01 ## [2,] 1 -4.12 0.01 ## [3,] 2 -3.82 0.01 ## [4,] 3 -3.36 0.01 ## [5,] 4 -3.00 0.01 ## Type 2: with drift no trend ## lag ADF p.value ## [1,] 0 -8.47 0.01 ## [2,] 1 -6.03 0.01 ## [3,] 2 -5.88 0.01 ## [4,] 3 -5.41 0.01 ## [5,] 4 -5.17 0.01 ## Type 3: with drift and trend ## lag ADF p.value ## [1,] 0 -8.46 0.01 ## [2,] 1 -6.02 0.01 ## [3,] 2 -5.86 0.01 ## [4,] 3 -5.38 0.01 ## [5,] 4 -5.14 0.01 ## ---- ## Note: in fact, p.value = 0.01 means p.value <= 0.01
3.3 短期利率序列单整性检验
1.平稳性检验
运行程序：
adf.test(SF) #短期利率序列平稳性检验
运行结果：
## Augmented Dickey-Fuller Test ## alternative: stationary ## ## Type 1: no drift no trend ## lag ADF p.value ## [1,] 0 -0.512 0.496 ## [2,] 1 -0.687 0.433 ## [3,] 2 -0.373 0.536 ## [4,] 3 -0.579 0.471 ## [5,] 4 -0.528 0.490 ## Type 2: with drift no trend ## lag ADF p.value ## [1,] 0 -1.96 0.346 ## [2,] 1 -2.27 0.225 ## [3,] 2 -1.86 0.384 ## [4,] 3 -2.13 0.278 ## [5,] 4 -2.03 0.316 ## Type 3: with drift and trend ## lag ADF p.value ## [1,] 0 -2.24 0.474 ## [2,] 1 -2.85 0.220 ## [3,] 2 -1.91 0.610 ## [4,] 3 -2.66 0.298 ## [5,] 4 -2.55 0.343 ## ---- ## Note: in fact, p.value = 0.01 means p.value <= 0.01
2.一阶差分后平稳性检验
运行程序：
adf.test(diff(SF)) #一阶差分短期利率序列平稳性检验
运行结果：
## Augmented Dickey-Fuller Test ## alternative: stationary ## ## Type 1: no drift no trend ## lag ADF p.value ## [1,] 0 -9.79 0.01 ## [2,] 1 -10.03 0.01 ## [3,] 2 -6.07 0.01 ## [4,] 3 -5.61 0.01 ## [5,] 4 -4.15 0.01 ## Type 2: with drift no trend ## lag ADF p.value ## [1,] 0 -9.77 0.01 ## [2,] 1 -10.02 0.01 ## [3,] 2 -6.07 0.01 ## [4,] 3 -5.61 0.01 ## [5,] 4 -4.15 0.01 ## Type 3: with drift and trend ## lag ADF p.value ## [1,] 0 -9.78 0.01 ## [2,] 1 -10.06 0.01 ## [3,] 2 -6.10 0.01 ## [4,] 3 -5.64 0.01 ## [5,] 4 -4.18 0.01 ## ---- ## Note: in fact, p.value = 0.01 means p.value <= 0.01
3.4 长期利率序列单整性检验
1.平稳性检验
运行程序：
adf.test(LF) #长期利率序列平稳性检验
运行结果：
## Augmented Dickey-Fuller Test ## alternative: stationary ## ## Type 1: no drift no trend ## lag ADF p.value ## [1,] 0 0.783 0.867 ## [2,] 1 0.451 0.772 ## [3,] 2 0.454 0.773 ## [4,] 3 0.297 0.728 ## [5,] 4 0.324 0.736 ## Type 2: with drift no trend ## lag ADF p.value ## [1,] 0 -1.08 0.667 ## [2,] 1 -1.29 0.594 ## [3,] 2 -1.30 0.591 ## [4,] 3 -1.40 0.557 ## [5,] 4 -1.32 0.584 ## Type 3: with drift and trend ## lag ADF p.value ## [1,] 0 -1.72 0.691 ## [2,] 1 -2.32 0.438 ## [3,] 2 -2.34 0.431 ## [4,] 3 -3.01 0.156 ## [5,] 4 -3.00 0.161 ## ---- ## Note: in fact, p.value = 0.01 means p.value <= 0.01
2.一阶差分后平稳性检验
运行程序：
adf.test(diff(LF)) #长阶差分短期利率序列平稳性检验
运行结果：
## Augmented Dickey-Fuller Test ## alternative: stationary ## ## Type 1: no drift no trend ## lag ADF p.value ## [1,] 0 -8.98 0.01 ## [2,] 1 -7.02 0.01 ## [3,] 2 -4.96 0.01 ## [4,] 3 -4.68 0.01 ## [5,] 4 -5.10 0.01 ## Type 2: with drift no trend ## lag ADF p.value ## [1,] 0 -9.05 0.01 ## [2,] 1 -7.11 0.01 ## [3,] 2 -5.05 0.01 ## [4,] 3 -4.77 0.01 ## [5,] 4 -5.22 0.01 ## Type 3: with drift and trend ## lag ADF p.value ## [1,] 0 -9.02 0.01 ## [2,] 1 -7.09 0.01 ## [3,] 2 -5.03 0.01 ## [4,] 3 -4.74 0.01 ## [5,] 4 -5.20 0.01 ## ---- ## Note: in fact, p.value = 0.01 means p.value <= 0.01
根据各变量的单位根检验，发现，这四个变量原始时序数据都存在单位根，一阶差分后平稳，表明这四个序列均为一阶单整序列。
4 各序列ARIMA模型拟合
4.1 log（M1）的ARIMA模型拟合
4.1.1 模型定阶
运行程序：
par(mfrow=c(1,2)) #创建1×2画布 M1<-na.omit(M1) #去除空值 acf(diff(M1),sub="图5 log（M）自相关图") #自相关图 pacf(diff(M1),sub="图6 log（M）偏自相关图") #偏自相关图
运行结果：
结合log（M）序列的ACF、Pacf特征,对其建立ARIMA（1，1，0）模型。
4.1.2 模型拟合
运行程序：
fit1<-arima(M1,order = c(1,1,0)) #拟合ARIMA（1，1，0）模型 fit1
运行结果：
## ## Call: ## arima(x = M1, order = c(1, 1, 0)) ## ## Coefficients: ## ar1 ## 0.5811 ## s.e. 0.0699 ## ## sigma^2 estimated as 0.0001071: log likelihood = 422.15, aic = -840.29
4.1.3 残差检验
运行程序：
ts.diag(fit1) #残差检验
运行结果：
根据结果显示，拟合效果较好。且残差序列为白噪声序列。该ARIMA模型为：∇ln(Mt)=1−0.5811B1ϵt
4.2 log（GNP）的ARIMA模型拟合
4.2.1 模型定阶
运行程序：
par(mfrow=c(1,2)) #创建1×2画布 acf(diff(GNP),sub="图8 log（GNP）自相关图") #自相关图 pacf(diff(GNP),sub="图9 log（GNP）偏自相关图") #偏自相关图
运行结果：
结合log（M）序列的ACF、Pacf特征,对其建立ARIMA（1，1，2）模型。
4.2.2 模型拟合
运行程序：
fit2<-arima(GNP,order = c(1,1,2)) #拟合ARIMA（1，1，2）模型 fit2
运行结果：
## ## Call: ## arima(x = GNP, order = c(1, 1, 2)) ## ## Coefficients: ## ar1 ma1 ma2 ## 0.9999 -0.7650 -0.2249 ## s.e. 0.0004 0.0722 0.0706 ## ## sigma^2 estimated as 9.719e-05: log likelihood = 430.54, aic = -853.09
4.2.3 残差检验
运行程序：
ts.diag(fit2) #残差检验
运行结果：
根据结果显示，拟合效果较好。且残差序列为白噪声序列。该ARIMA模型为：
∇ln(GNPt)=1−0.9999B1−0.765B−0.2249B2ϵt
4.3 短期利率序列的ARIMA模型拟合
4.3.1 模型定阶
运行程序：
par(mfrow=c(1,2)) #创建1×2画布 acf(diff(SF),sub="图11 短期利率自相关图") #自相关图 pacf(diff(SF),sub="图12 短期利率偏自相关图 #偏自相关图
运行结果：
结合log（M）序列的ACF、Pacf特征,对其建立ARIMA（7，1，0）模型。
4.4.2 模型拟合
运行程序：
fit3<-arima(SF,order = c(7,1,0)) #拟合ARIMA（7，1，0）模型 fit3
运行结果：
## ## Call: ## arima(x = SF, order = c(7, 1, 0)) ## ## Coefficients: ## ar1 ar2 ar3 ar4 ar5 ar6 ar7 ## 0.2788 -0.3515 0.2863 -0.0879 0.1233 -0.0918 -0.2366 ## s.e. 0.0832 0.0860 0.0903 0.0930 0.0896 0.0847 0.0819 ## ## sigma^2 estimated as 5.599e-05: log likelihood = 468.76, aic = -921.53
4.4.4 疏系数模型拟合
运行程序：
fit31<-arima(SF,order = c(7,1,0),transform.pars=F,fixed=c(NA,NA,NA,0,0,0,NA)) #拟合ARIMA（7，1，0）疏系数模型 fit31
运行结果：
## ## Call: ## arima(x = SF, order = c(7, 1, 0), transform.pars = F, fixed = c(NA, NA, NA, ## 0, 0, 0, NA)) ## ## Coefficients: ## ar1 ar2 ar3 ar4 ar5 ar6 ar7 ## 0.2477 -0.3080 0.2223 0 0 0 -0.2907 ## s.e. 0.0796 0.0776 0.0792 0 0 0 0.0742 ## ## sigma^2 estimated as 5.699e-05: log likelihood = 467.61, aic = -925.22
4.3.4 残差检验
运行程序：
ts.diag(fit31) #残差检验
运行结果：
根据结果显示，拟合效果较好。且残差序列为白噪声序列。该ARIMA模型为：
∇SFt=1−0.2477B+0.3080B2−0.2223B3+0.2907B71ϵt
4.4 长期利率序列的ARIMA模型拟合
4.4.1 模型定阶
运行程序：
library(forecast) par(mfrow=c(1,2)) #创建1×2画布 acf(diff(LF),sub="图14 长期利率自相关图") #自相关图 pacf(diff(LF),sub="图15 长期利率自相关图") #偏自相关图
运行结果：
结合log（M）序列的ACF、Pacf特征,对其建立ARIMA（14，1，0）模型。
4.4.2 模型拟合
运行程序：
fit4<-arima(LF,order=c(14,1,0)) #拟合ARIMA（14，1，0）模型 fit4
运行结果：
## ## Call: ## arima(x = LF, order = c(14, 1, 0)) ## ## Coefficients: ## ar1 ar2 ar3 ar4 ar5 ar6 ar7 ar8 ## 0.2619 -0.0558 0.1885 -0.0039 -0.1647 -0.0194 -0.0940 0.0245 ## s.e. 0.0815 0.0836 0.0843 0.0853 0.0848 0.0852 0.0846 0.0871 ## ar9 ar10 ar11 ar12 ar13 ar14 ## 0.0934 0.0383 -0.1036 0.1030 -0.1081 0.2992 ## s.e. 0.0871 0.0859 0.0859 0.0849 0.0873 0.0854 ## ## sigma^2 estimated as 1.392e-05: log likelihood = 562.28, aic = -1094.56
4.4.3 疏系数模型拟合
运行程序：
fit41<-arima(LF,order = c(14,1,0),transform.pars=F,fixed=c(NA,0,0,0,NA,0,0,0,0,0,0,0,0,NA)) #拟合ARIMA（14，1，0）疏系数模型 fit41
运行结果：
## ## Call: ## arima(x = LF, order = c(14, 1, 0), transform.pars = F, fixed = c(NA, 0, 0, 0, ## NA, 0, 0, 0, 0, 0, 0, 0, 0, NA)) ## ## Coefficients: ## ar1 ar2 ar3 ar4 ar5 ar6 ar7 ar8 ar9 ar10 ar11 ar12 ## 0.2392 0 0 0 -0.1675 0 0 0 0 0 0 0 ## s.e. 0.0785 0 0 0 0.0794 0 0 0 0 0 0 0 ## ar13 ar14 ## 0 0.2566 ## s.e. 0 0.0835 ## ## sigma^2 estimated as 1.516e-05: log likelihood = 556.84, aic = -1105.69
4.4.4 残差检验
运行程序：
ts.diag(fit41) #残差检验
运行结果：
根据结果显示，拟合效果较好。且残差序列为白噪声序列。该ARIMA模型为：
∇LFt=1−0.2392B+0.1672B5−0.2566B141ϵt
