VAR模型分析聯合內生變數的動態關係

一、實驗介紹

1.1 實驗內容

VAR模型是向量自迴歸模型的簡稱，是基於資料的統計性質建立的一種常用的計量經濟模型，它把系統中每一個內生變數作為系統中所有內生變數的滯後值的函式來構造模型，從而將單變數自迴歸模型推廣到由多元時間序列變數組成的“向量”自迴歸模型。本實驗運用 R 語言來建立兩變數的向量自迴歸模型，首先是檢驗兩變數序列的平穩性，然後進行協整檢驗，確定待擬合模型的滯後階數，再擬合VAR模型和做脈衝響應分析，最終對擬合的VAR模型進行分析於預測。通過本實驗學會用VAR模型處理多個相關經濟指標的分析與預測。

1.2 實驗知識點

• 平穩性檢驗
• 協整檢驗
• 滯後階數的確定
• VAR 模型的擬合
• 脈衝響應分析
• VAR 模型的預測

1.3 實驗環境

• R version 3.4.1
• Xfce 終端

1.4 適合人群

本課程難度為較難，屬於中級級別課程，適合具有 R 語言基礎並且有一定經管背景的使用者，熟悉 R 語言基礎知識並加深鞏固，學會用 R 語言建立計量經濟模型分析經濟資料。

1.5 實驗準備

開啟Xfce 終端，下載實驗所需的資料，並啟動 R ：

$ wget http://labfile.oss.aliyuncs.com/courses/910/monthdata.csv
$ sudo R

二、實驗原理

以下實驗原理介紹來自 MBA 智庫百科：向量自迴歸模型

向量自迴歸模型簡稱VAR模型，是一種常用的計量經濟模型，1980年由克里斯托弗·西姆斯（Christopher Sims）提出。VAR模型是用模型中所有當期變數對所有變數的若干滯後變數進行迴歸。VAR模型用來估計聯合內生變數的動態關係，而不帶有任何事先約束條件。它是AR模型的推廣，此模型目前已得到廣泛應用。

向量自迴歸(VAR)是基於資料的統計性質建立模型，VAR模型把系統中每一個內生變數作為系統中所有內生變數的滯後值的函式來構造模型，從而將單變數自迴歸模型推廣到由多元時間序列變數組成的“向量”自迴歸模型。VAR模型是處理多個相關經濟指標的分析與預測最容易操作的模型之一，並且在一定的條件下，多元MA和ARMA模型也可轉化成VAR模型，因此近年來VAR模型受到越來越多的經濟工作者的重視。

三、實驗步驟

本實驗利用 2005 年 8 月匯率改革後消費者價格指數（CPI）和人民幣名義有效匯率（NEER）的月度資料建立兩變數 VAR 模型來分析人民幣匯率變動與CPI之間的動態關係。

3.1 實驗準備

> data<-read.csv("monthdata.csv"
 
)  #讀取資料
> CPI<-data[,1]    #將 data 資料的第一列賦值給 CPI 資料框
> NEER<-data[,2]    
> length(CPI)  #檢視 CPI 變數的資料長度
> length(NEER)

可以看到資料長度為 137，即包含了 137 個月度 CPI 和 NEER 資料。

下面我們繪製兩個變數的時序圖，觀察消費者價格指數和名義匯率的波動情況。

> CPI.ts<-ts(CPI,start=c(2005,8),end=c(2016,12),freq=12)   #將資料轉化為時間序列
> NEER.ts<-ts(NEER,start=c(2005,8),end=c(2016,12),freq=12)
> par(mfrow=c(2,1))
> plot(CPI.ts,type="l",xlab="Date", ylab="CPI")     
> plot(NEER.ts,type="l",xlab="Date", ylab="NEER")

可以看到消費者價格指數（CPI）的波動幅度較大，而名義有效匯率波動較小但明顯呈現一個波動上升的趨勢，兩個變數序列大致是不平穩的。

3.2 平穩性檢驗

在擬合 VAR 模型之前，需要對變數進行平穩性檢驗，如果要擬合的內生變數都是平穩的或者同階單整的才可進行 VAR 模型的擬合。但在平穩性檢驗之前需要先項兩變數取對數，以消除時間序列的異方差的影響。

> lncpi<-log(CPI)     
> lnneer<-log(NEER)

進行平穩性檢驗本實驗選擇 urca 包中的 ur.df 函式進行單位根檢驗，若存在單位根則序列不平穩；反之，不存在單位根則序列平穩。

> install.packages("urca")    #映象選擇 China (Lanzhou) [https]
> library(urca)
> urt.lncpi<-ur.df(lncpi,type='trend',selectlags='AIC')
> urt.lnneer<-ur.df(lnneer,type='trend',selectlags='AIC')
> summary(urt.lncpi)    

  ############################################### 
  # Augmented Dickey-Fuller Test Unit Root Test # 
  ############################################### 

  Test regression trend 


  Call:
  lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)

  Residuals:
        Min        1Q    Median        3Q       Max 
  -0.028676 -0.003311  0.000381  0.003035  0.043118 

  Coefficients:
                Estimate Std. Error t value Pr(>|t|)  
  (Intercept)  3.047e-01  1.477e-01   2.063   0.0411 *
  z.lag.1     -6.555e-02  3.185e-02  -2.058   0.0416 *
  tt          -1.443e-05  1.678e-05  -0.860   0.3912  
  z.diff.lag  -1.178e-01  8.627e-02  -1.365   0.1745  
  ---
  Signif. codes:  0 \u2018***\u2019 0.001 \u2018**\u2019 0.01 \u2018*\u2019 0.05 \u2018.\u2019 0.1 \u2018 \u2019 1

  Residual standard error: 0.007505 on 131 degrees of freedom
  Multiple R-squared:  0.05455,    Adjusted R-squared:  0.0329 
  F-statistic: 2.519 on 3 and 131 DF,  p-value: 0.06084


  Value of test-statistic is: -2.0578 1.5223 2.2719 

  Critical values for test statistics: 
        1pct  5pct 10pct
  tau3 -3.99 -3.43 -3.13
  phi2  6.22  4.75  4.07
  phi3  8.43  6.49  5.47

> summary(urt.lnneer)   

  ############################################### 
  # Augmented Dickey-Fuller Test Unit Root Test # 
  ############################################### 

  Test regression trend 


  Call:
  lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)

  Residuals:
        Min     1Q    Median        3Q       Max 
  -0.038574 -0.007164 -0.000990  0.009182  0.033482 

  Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
  (Intercept)  0.4429099  0.1380821   3.208  0.00168 ** 
  z.lag.1     -0.0995380  0.0312008  -3.190  0.00178 ** 
  tt           0.0003157  0.0001058   2.983  0.00341 ** 
  z.diff.lag   0.4099249  0.0810914   5.055 1.42e-06 ***
  ---
  Signif. codes:  0 \u2018***\u2019 0.001 \u2018**\u2019 0.01 \u2018*\u2019 0.05 \u2018.\u2019 0.1 \u2018 \u2019 1

  Residual standard error: 0.01267 on 131 degrees of freedom
  Multiple R-squared:  0.1877,    Adjusted R-squared:  0.1691 
  F-statistic: 10.09 on 3 and 131 DF,  p-value: 5.01e-06


  Value of test-statistic is: -3.1902 4.2783 5.1509 

  Critical values for test statistics: 
        1pct  5pct 10pct
  tau3 -3.99 -3.43 -3.13
  phi2  6.22  4.75  4.07
  phi3  8.43  6.49  5.47

判斷是否平穩主要看詳細擬合結果的最後兩個部分，即：

• 1.Value of test-statistic is:
• 2.Critical values for test statistics:

1 是檢驗統計量的值，2 是對應的顯著性水平下檢驗統計量的臨界值。

單位根檢驗的原假設是序列存在單位根。

對於 lncpi 和 lnneer 的單位根檢驗結果:

• lncpi的檢驗統計量的值 -2.0578 在 1%、5%、10% 的顯著性水平下都大於臨界值(-3.99 -3.43 -3.13)，則不能拒絕原價設而接受存在單位根的假設，說明 lncpi 序列是不平穩的。（注意，統計量的值和臨界值為負，統計量的值大於臨界值是接受原假設；若統計量的值和臨界值為正值，統計量的值大於臨界值是拒絕原假設。）
• 同理，lnneer 的檢驗統計量的值 -3.1902 在 1% 、5% 的顯著性水平下都大於對應的臨界值(-3.99 -3.43 -3.13)，不能拒絕原假設，lnneer 序列也是不平穩的，存在單位根。

由於兩個變數都存在單位根，對序列差分後檢驗序列的平穩性。

> dlncpi<-diff(lncpi)   #取一階差分
> dlnneer<-diff(lnneer)
> urt.dlncpi<-ur.df(dlncpi,type='trend',selectlags='AIC')
> urt.dlnneer<-ur.df(dlnneer,type='trend',selectlags='AIC')
> summary(urt.dlncpi)  
> summary(urt.dlnneer)

由於檢驗的詳細結果太長，擷取需要的部分一階差分平穩性檢驗結果：

• lncpi 的一階差分平穩性檢驗

  

  由於 dlncpi 的檢驗統計量的值為 -8.1895，在 1%、5%、10% 的顯著性水平上都小於對應的臨界值，因此拒絕存在單位根的原假設，即 lncpi 是一階差分平穩的。

• lnneer 的一階差分平穩性檢驗

  

  由於 dlnneer 的檢驗統計量的值為 -7.0409，在 1%、5%、10% 的顯著性水平上都小於對應的臨界值，因此拒絕存在單位根的原假設，即 lncpi 是一階差分平穩的。

我們再來看一階差分後的時間序列圖：

> par(mfrow=c(2,1))
> plot(dlncpi,type="l",xlab="Date", ylab="diff.CPI")     
> plot(dlnneer,type="l",xlab="Date", ylab="diff.NEER")

可以看到差分後的兩個變數序列沒有明顯的波動聚集或則上升下降的趨勢，比較平穩。

兩個內生變數都是一階差分平穩的，即都是一階單整的，是不平穩的時間序列，不能做 Granger 因果檢驗，只能做協整檢驗。

3.3 協整檢驗

協整檢驗主要針對非平穩的單個序列，但它們的線性組合可能是平穩的。幾個變數之間可能存在的一種長期均衡關係進行檢驗，表現為存在某個協整方程。

由於所有變數都是一階單整的，是非平穩時間序列，因此各變數之間可能存在協整關係，如果要對所選擇的內生變數進行VAR模型的構建，需要進行協整檢驗，以判斷各個變數之間是否存在長期穩定的協整關係，處理各變數之間的是否存在偽迴歸問題。在本實驗我們運用E-G兩步法協整檢驗，需要用到 lmtest 包中的 dwtest 函式。

• E-G兩步法：第一步:迴歸方程的估計

> fit<-lm(lncpi~lnneer)  
> summary(fit)

由於擬合的模型變數和截距項的 t 檢驗在 5% 的顯著性水平上都是顯著的，且模型的 F 檢驗在5%的顯著性水平上也是顯著的，因此擬合的線性模型是合理的。

協整迴歸方程為：lncpi=4.80296-0.03666 lnneer+εt

> install.packages("lmtest")
> library(zoo)
> library(lmtest)
> dwtest(fit)    #檢驗序列的自相關性

由於 Durbin-Watson test 檢驗的 p-value 的值幾乎為 0，該數值小於 0.05 ，說明在 5% 的顯著性水平上殘差序列不獨立，具有自相關性。

• 檢驗殘差序列的平穩性：

  > error<-residuals(fit)  #提取殘差序列
  > urt.residuals<-ur.df(error,type='none',selectlags='AIC')
  > summary(urt.residuals)
  

  

  

  根據殘差序列的平穩性檢驗結果，在 5% 的顯著性水平上拒絕殘差序列存在單位根的原假設，即殘差序列是平穩的，說明 CPI 和名義有效匯率兩個序列之間存在協整關係，意味著我國的消費者價格指數和名義有效匯率之間具有長期均衡關係，增長或者減少具有協同效應。

  。

• E-G兩步法：第二步：誤差修正模型的建立

> error.lag<-error[-c(137,138)]       
> ecm.fit<-lm(dlncpi~error.lag+dlnneer)   #擬合誤差修正模型
> summary(ecm.fit)  
> dwtest(ecm.fit)

協整迴歸方程：

誤差修正項的係數為負，符合誤差修正機制，反映了上一期偏離長期均衡的數量將在下一期得到反向修正，這也符合之前證明的協整關係。

3.4 滯後階數的確定

在擬合 VAR 模型之前還需要確定擬合幾階 VAR 模型，也就是確定滯後階數。確定滯後階數需要用到 vars 包中的 VARselect 函式。

> install.packages("vars")
> library(MASS)
> library(sandwich)
> library(strucchange)
> library(vars)
> data.new<-data.frame(lncpi,lnneer)  #合併資料
> VARselect(data.new,lag.max=10,type="const")  #在10以內選擇最優滯後階數

根據結果，不同的資訊準則有不同的滯後階數，選擇 2 階或者 4 階都是可以的，一般來說選擇在相同條件下更加簡潔的模型，因此選擇 2 階滯後。

3.5 VAR模型的擬合和預測

3.5.1.VAR 模型的擬合

在確定好最優滯後階數以後我們就可以擬合模型：

> var<-VAR(data.new,lag.max=2,ic="AIC")
> summary(var)

  VAR Estimation Results:   #VAR(2) 模型擬合結果
  ========================= 
  Endogenous variables: lncpi, lnneer 
  Deterministic variables: const 
  Sample size: 135 
  Log Likelihood: 873.227 
  Roots of the characteristic polynomial:
  0.9621 0.9621 0.3829 0.1847
  Call:
  VAR(y = data.new, lag.max = 2, ic = "AIC")

  Estimation results for equation lncpi:     #模型估計結果一
  ====================================== 
  lncpi = lncpi.l1 + lnneer.l1 + lncpi.l2 + lnneer.l2 + const 

            Estimate Std. Error t value Pr(>|t|)    
  lncpi.l1   0.80917    0.08578   9.433   <2e-16 ***
  lnneer.l1 -0.09248    0.04705  -1.965   0.0515 .  
  lncpi.l2   0.12907    0.08483   1.521   0.1306    
  lnneer.l2  0.08381    0.04712   1.779   0.0777 .  
  const      0.32677    0.15855   2.061   0.0413 *  
  ---
  Signif. codes:  0 \u2018***\u2019 0.001 \u2018**\u2019 0.01 \u2018*\u2019 0.05 \u2018.\u2019 0.1 \u2018 \u2019 1


  Residual standard error: 0.007369 on 130 degrees of freedom
  Multiple R-Squared: 0.8791,    Adjusted R-squared: 0.8754 
  F-statistic: 236.3 on 4 and 130 DF,  p-value: < 2.2e-16 


  Estimation results for equation lnneer:        #模型估計結果二
  ======================================= 
  lnneer = lncpi.l1 + lnneer.l1 + lncpi.l2 + lnneer.l2 + const 

            Estimate Std. Error t value Pr(>|t|)    
  lncpi.l1  -0.18430    0.14954  -1.232 0.220030    
  lnneer.l1  1.31275    0.08203  16.003  < 2e-16 ***
  lncpi.l2   0.28850    0.14789   1.951 0.053238 .  
  lnneer.l2 -0.31998    0.08215  -3.895 0.000156 ***
  const     -0.44726    0.27641  -1.618 0.108065    
  ---
  Signif. codes:  0 \u2018***\u2019 0.001 \u2018**\u2019 0.01 \u2018*\u2019 0.05 \u2018.\u2019 0.1 \u2018 \u2019 1


  Residual standard error: 0.01285 on 130 degrees of freedom
  Multiple R-Squared: 0.9908,    Adjusted R-squared: 0.9905 
  F-statistic:  3504 on 4 and 130 DF,  p-value: < 2.2e-16 


  Covariance matrix of residuals:       #殘差的協方差矩陣
             lncpi    lnneer
  lncpi  5.430e-05 7.761e-06
  lnneer 7.761e-06 1.650e-04

  Correlation matrix of residuals:      #殘差的相關陣
           lncpi  lnneer
  lncpi  1.00000 0.08198
  lnneer 0.08198 1.00000

根據模型擬合結果，所得的 VAR(2) 的模型方程（保留兩位小數）為：

也可以用 coef 函式來檢視模型估計的簡明結果。

> coef(var)   #檢視模型估計的簡明結果

會擬合出一致的 VAR(2) 方程（係數保留兩位小數）。

> plot(var)       #畫出每個變數的時序圖、殘差圖、ACF 圖、PACF 圖
按<Return>鍵來看下一個圖      #也就是<Enter>鍵

• 變數 lncpi 的時序圖、殘差圖、ACF 圖、PACF 圖
• 變數 lnneer 的時序圖、殘差圖、ACF 圖、PACF 圖

3.5.2 脈衝響應分析

在檢視完擬合結果的圖形之後，我們來繪製擬合結果的脈衝響應圖，需要用到 var 包中的 irf 函式。然後通過 plot 函式畫出圖形。

> var.irf<-irf(var)  
> plot(var.irf)       
按<Return>鍵來看下一個圖：

根據圖形 lncpi 自身以及 lnneer 的波動對 lncpi 有正向的衝擊。 lncpi 對自身的影響沒有滯後期，並且自身波動的影響隨著時期的增加會越來越小。lnneer 波動對 lncpi 的影響在第一期以前是逐漸減少的，但在第一期減少為 0 之後隨著時間的增加影響越來越大。

lncpi 波動對 lnneer 有一個負向的衝擊，並且隨著時期數的增加負向的影響會越來越大。lnneer 波動對自身有一個正向的衝擊，這個衝擊從開始先增加，在第二期達到最大值以後又逐漸減少。

3.5.3 VAR(2)模型的預測

模型預測通常需要用到 predict 函式，具體用法如下：

> var.predict<-predict(var,n.ahead=10,ci=0.95)
> var.predict

這樣我們就得到了兩變數（lncpi、lnneer）的 VAR(2)模型的滯後 10 期預測結果：

四、實驗總結

通過本次實驗深入瞭解向量自迴歸模型，學會運用 R 語言來進行 VAR模型建模。學習如何檢驗一個序列的平穩性，以及如何運用協整檢驗來分析非平穩時間序列的關係。學會運用 var 包中的包括 VARselect() 函式、VAR() 函式、coef() 函式以及 irf() 函式各種函式來進行 VAR 模型的建模，選擇最優滯後階數，擬合 VAR 模型，對模型結果進行脈衝響應分析以及模型的預測。希望同過實驗能夠獨立的運用 R 語言進行 VAR 模型建模於分析，鞏固 R 的同時加深自己的金融知識。