1 多元时间序列分析

1.1 平稳多元序列建模:动态回归模型ARIMAX

# 这是以中心化后的天然气输入速率为输入序列,建立CO2的输出百分浓度模型
# 读入数据,并绘制输出序列时序图
a<-read.table("file24.csv",sep=",",header = T)
y<-ts(a$output)
plot(y)

# 考察输出序列自相关图和偏自相关图
acf(y)

pacf(y)

# 对输出序列拟合模型AR(4)
y.fit<-arima(y,order=c(4,0,0))
y.fit
## 
## Call:
## arima(x = y, order = c(4, 0, 0))
## 
## Coefficients:
##          ar1      ar2      ar3     ar4  intercept
##       2.0976  -1.3243  -0.0071  0.2124    53.6880
## s.e.  0.0565   0.1344   0.1345  0.0566     0.8679
## 
## sigma^2 estimated as 0.1105:  log likelihood = -97.28,  aic = 206.56
# 残差白噪声检验
for (i in 1:2) print(Box.test(y.fit$residual,lag=6*i))
## 
##  Box-Pierce test
## 
## data:  y.fit$residual
## X-squared = 8.3004, df = 6, p-value = 0.2169
## 
## 
##  Box-Pierce test
## 
## data:  y.fit$residual
## X-squared = 12.811, df = 12, p-value = 0.3829

x<-ts(a$input)
ccf(y,x)

R语言中的 TSA 包中的 ‘arimax’ 函数可以拟合 ARIMAX 模型,

library(TSA)
## 
## Attaching package: 'TSA'
## The following objects are masked from 'package:stats':
## 
##     acf, arima
## The following object is masked from 'package:utils':
## 
##     tar
y.fit2<-arimax(y,order=c(4,0,0),xtransf = x,transfer = list(c(2,1)))
y.fit2
## 
## Call:
## arimax(x = y, order = c(4, 0, 0), xtransf = x, transfer = list(c(2, 1)))
## 
## Coefficients:
##          ar1      ar2      ar3     ar4  intercept  T1-AR1   T1-AR2  T1-MA0
##       1.6469  -0.8310  -0.0639  0.2095    53.5836  1.5508  -0.6667  0.3122
## s.e.  0.0585   0.1119   0.1129  0.0595     0.4030  0.0275   0.0266  0.0621
##        T1-MA1
##       -0.5943
## s.e.   0.0677
## 
## sigma^2 estimated as 0.07143:  log likelihood = -31.25,  aic = 80.49

1.2 虚假回归 / 伪回归

1.3 单位根检验

1.3.1 非平稳序列类型

下面拟合 \(\Phi_0=0.1\) 带漂移项的差分平稳序列 \(x_t=\Phi_0+x_{t-1}+\epsilon_t\) ,并绘制时序图和线性拟合效果图

e<-rnorm(1000)
x<-filter(0.1+e,filter = 1,method = "recursive")
plot(x)
t<-c(1:1000)
abline(lm(x~t),col=2)

lm(x~t)
## 
## Call:
## lm(formula = x ~ t)
## 
## Coefficients:
## (Intercept)            t  
##     2.74115      0.02607
#差分运算和线性拟合后残差序列比较
r1<-diff(x)
x.fit<-lm(x~t)
r2<-ts(x.fit$residual)
c1<-min(r1,r2)
c2<-max(r1,r2)
plot(r1,ylim=c(c1,c2))
lines(r2,col=2)

对于带漂移项的差分平稳序列(差分后序列平稳),尽管时序图显示该序列具有显著的线性趋势,但是用差分运算提取趋势信息要比用线性拟合模型提取线性趋势更加清晰。

下面拟合 \(\Phi_0=0, \beta=0.1\) 趋势平稳序列 \(x_t=\Phi_0+\beta t+\epsilon_t\) ,并绘制时序图

t<-c(1:1000)
e<-rnorm(1000,0,10)
x<-0.1*t+e
x<-ts(x)
plot(x)

#考察差分后差分序列的波动范围
plot(diff(x))
abline(h=c(1.98*sd(diff(x)),-1.98*sd(diff(x))),col=2)

#考察线性拟合后残差序列的波动范围
x.fit<-lm(x~t)
r<-ts(x.fit$residual)
plot(r)
abline(h=c(1.98*sd(r),-1.98*sd(r)),col=2)

对于趋势平稳序列,尽管差分运算与趋势拟合均能实现序列平稳,但是趋势拟合的残差波动范围更窄,拟合程度更高。

1.3.2 DF检验:fUnitRoots 程序包 adfTest

下面是1978-2002年中国农村居民家庭人均纯收入对数序列 \(ln(x_t)\) 与消费支出对数序列 \(ln(y_t)\) 进行DF检验

#读入数据,并绘制两序列时序图
b<-read.table("file25.csv",sep=",",header = T)
x<-ts(b$lnx,start = 1978)
y<-ts(b$lny,start = 1978)
c1<-min(x,y)
c2<-max(x,y)
plot(x,ylim=c(c1,c2))
lines(y,lty=2,col=2)

#对人均纯收入对数序列进行DF检验
library(fBasics)
## Loading required package: timeDate
## 
## Attaching package: 'timeDate'
## The following objects are masked from 'package:TSA':
## 
##     kurtosis, skewness
## Loading required package: timeSeries
library(fUnitRoots)
adfTest(x,lag=1,type="nc")
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: 1.0859
##   P VALUE:
##     0.9202 
## 
## Description:
##  Sun Dec 19 18:09:04 2021 by user: 83930
adfTest(x,lag=1,type="c")
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -1.1606
##   P VALUE:
##     0.6201 
## 
## Description:
##  Sun Dec 19 18:09:04 2021 by user: 83930
adfTest(x,lag=1,type="ct")
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -3.1381
##   P VALUE:
##     0.1388 
## 
## Description:
##  Sun Dec 19 18:09:04 2021 by user: 83930
#对人均生活消费支出对数序列进行DF检验
adfTest(y,lag=1,type="nc")
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: 1.0806
##   P VALUE:
##     0.9196 
## 
## Description:
##  Sun Dec 19 18:09:04 2021 by user: 83930
adfTest(y,lag=1,type="c")
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -1.4086
##   P VALUE:
##     0.5324 
## 
## Description:
##  Sun Dec 19 18:09:04 2021 by user: 83930
adfTest(y,lag=1,type="ct")
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -3.3125
##   P VALUE:
##     0.08993 
## 
## Description:
##  Sun Dec 19 18:09:04 2021 by user: 83930

显著性水平取 0.05 时,可认为 中国农村居民家庭人均纯收入对数序列 与 消费支出对数序列 均为非平稳序列,这与时序图显示出来的性质一致。

1.3.3 ADF检验

由于DF检验只适用于 \(AR(1)\) 的平稳性检验,但实际上绝大多数时间序列都不会是一个简单的 \(AR(1)\) 过程。为适应于 \(AR(p)\) 过程的平稳性检验,提出 ADF 检验。

下面是1978-2002年中国农村居民家庭人均纯收入对数序列差分后 \(\nabla ln(x_t)\) 与生活消费支出对数序列差分后 \(\nabla ln(y_t)\) 进行ADF检验

#差分运算
dx<-diff(x)
dy<-diff(y)

#人均收入对数差分后序列ADF检验
for (i in 1:3) print(adfTest(dx,lag=i,type='nc'))
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -1.3701
##   P VALUE:
##     0.173 
## 
## Description:
##  Sun Dec 19 18:09:04 2021 by user: 83930
## 
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 2
##   STATISTIC:
##     Dickey-Fuller: -1.3569
##   P VALUE:
##     0.1772 
## 
## Description:
##  Sun Dec 19 18:09:04 2021 by user: 83930
## 
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 3
##   STATISTIC:
##     Dickey-Fuller: -1.3367
##   P VALUE:
##     0.1836 
## 
## Description:
##  Sun Dec 19 18:09:04 2021 by user: 83930
for (i in 1:3) print(adfTest(dx,lag=i,type='c'))
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -2.2267
##   P VALUE:
##     0.2428 
## 
## Description:
##  Sun Dec 19 18:09:04 2021 by user: 83930
## 
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 2
##   STATISTIC:
##     Dickey-Fuller: -3.0963
##   P VALUE:
##     0.0427 
## 
## Description:
##  Sun Dec 19 18:09:04 2021 by user: 83930
## 
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 3
##   STATISTIC:
##     Dickey-Fuller: -1.8967
##   P VALUE:
##     0.3596 
## 
## Description:
##  Sun Dec 19 18:09:04 2021 by user: 83930
#人均支出对数差分后序列ADF检验
for (i in 1:3) print(adfTest(dy,lag=i,type='nc'))
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -1.7967
##   P VALUE:
##     0.0719 
## 
## Description:
##  Sun Dec 19 18:09:04 2021 by user: 83930
## 
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 2
##   STATISTIC:
##     Dickey-Fuller: -1.3048
##   P VALUE:
##     0.1937 
## 
## Description:
##  Sun Dec 19 18:09:04 2021 by user: 83930
## 
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 3
##   STATISTIC:
##     Dickey-Fuller: -1.0581
##   P VALUE:
##     0.272 
## 
## Description:
##  Sun Dec 19 18:09:04 2021 by user: 83930
for (i in 1:3) print(adfTest(dy,lag=i,type='c'))
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -3.3473
##   P VALUE:
##     0.02438 
## 
## Description:
##  Sun Dec 19 18:09:04 2021 by user: 83930
## 
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 2
##   STATISTIC:
##     Dickey-Fuller: -2.809
##   P VALUE:
##     0.0758 
## 
## Description:
##  Sun Dec 19 18:09:04 2021 by user: 83930
## 
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 3
##   STATISTIC:
##     Dickey-Fuller: -1.7169
##   P VALUE:
##     0.4232 
## 
## Description:
##  Sun Dec 19 18:09:04 2021 by user: 83930

结果显示,中国农村居民家庭人均纯收入对数序列差分后为具有常数均值的平稳序列,该平稳序列 2 阶自相关(p=0.0427<0.05),生活消费支出对数序列差分后同样为具有常数均值的平稳序列,该平稳序列 1 阶自相关(p=0.02438<0.05)

1.4 协整检验:Engle-Granger 检验

下面是1978-2002年中国农村居民家庭人均纯收入对数序列\(ln(x_t)\) 与生活消费支出对数序列 \(ln(y_t)\) 进行EG检验

#第一步:构造回归模型
y.fit<-lm(y~x)
summary(y.fit)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.083331 -0.043257 -0.003606  0.044968  0.083295 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.07361    0.07337   1.003    0.326    
## x            0.95729    0.01110  86.254   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05095 on 23 degrees of freedom
## Multiple R-squared:  0.9969, Adjusted R-squared:  0.9968 
## F-statistic:  7440 on 1 and 23 DF,  p-value: < 2.2e-16
#第二步:残差序列单位根检验
r<-ts(y.fit$residual)
for (i in 1:3) print(adfTest(r,lag=i,type='nc'))
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -1.8644
##   P VALUE:
##     0.06224 
## 
## Description:
##  Sun Dec 19 18:09:04 2021 by user: 83930
## 
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 2
##   STATISTIC:
##     Dickey-Fuller: -1.9873
##   P VALUE:
##     0.04699 
## 
## Description:
##  Sun Dec 19 18:09:04 2021 by user: 83930
## 
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 3
##   STATISTIC:
##     Dickey-Fuller: -1.581
##   P VALUE:
##     0.106 
## 
## Description:
##  Sun Dec 19 18:09:04 2021 by user: 83930

结果显示,尽管 中国农村居民家庭人均纯收入对数序列 与 生活消费支出对数序列 均为非平稳序列,但他们之间具有协整关系(回归残差序列属于无常数均值 2 阶自相关平稳序列,ADF检验 p=0.04699<0.05),收入与支出之间具有协整关系,故可建立如下回归模型拟合其长期均衡关系:

\[ \ln(y_t)=0.07361+0.95729 \ln(x_t) \]

1.5 误差修正模型:ECM

下面是1978-2002年中国农村居民家庭人均纯收入对数序列\(ln(x_t)\) 与生活消费支出对数序列 \(ln(y_t)\) 构造ECM模型

ECM<-y.fit$residual[1:24]
dify.fit<-lm(diff(y)~0+diff(x)+ECM)
summary(dify.fit)
## 
## Call:
## lm(formula = diff(y) ~ 0 + diff(x) + ECM)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.049691 -0.025353  0.003606  0.024944  0.049214 
## 
## Coefficients:
##         Estimate Std. Error t value Pr(>|t|)    
## diff(x)  0.95513    0.04472   21.36 3.37e-16 ***
## ECM     -0.17152    0.12796   -1.34    0.194    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0301 on 22 degrees of freedom
## Multiple R-squared:  0.9551, Adjusted R-squared:  0.951 
## F-statistic: 233.9 on 2 and 22 DF,  p-value: 1.501e-15

结果显示,构造 ECM 模型:

\[ \nabla y_t=0.95513 \nabla x_t-0.17152 ECM_{t-1} \]

收入当期波动对生活消费支出当期波动的影响很大,每增加 1 单位的对数收入,会增加 0.95513 单位的对数生活消费支出。

上期误差 ECM 对生活消费支出当期波动的调整幅度不大,单位调整比例为 -0.17152,且系数显著性检验显示该系数不显著非零。