1. 程式人生 > >Stata: VAR - 模擬、估計和推斷

Stata: VAR - 模擬、估計和推斷


作者:許夢潔 (編譯) (知乎 | 簡書 | 碼雲)

Stata連享會 「Stata 現場班報名中……」

Source: Ashish RajbhandariVector autoregression—simulation, estimation, and inference in Stata

文章目錄


VAR 是分析多維時間序列動態變化的利器,該模型設定為由一組時間序列組成的序列是其自己滯後項的函式。

1. 模擬

首先使用如下設定模擬雙變數 VAR(2) :

(1)

[ y 1
, t
y 2 , t
]
= μ + A 1 + [ y 1 , t 1 y 2 , t 1 ] + A 2 + [ y 1 , t 2 y 2 , t 2 ] + ϵ t
\begin{bmatrix} y_{ 1,t } \\ y_{ 2,t } \end{bmatrix}= \mu +A_{ 1 }+ \begin{bmatrix} y_{ 1,t-1 }\\ y_{ 2,t-1 } \end{bmatrix} +A_{ 2 }+ \begin{bmatrix} y_{ 1,t-2 }\\ y_{ 2,t-2 } \end{bmatrix} +\epsilon _{ t } \tag{1}

其中 y 1 , t y_{ 1,t } 是在時間 t t 的觀測變數, μ \mu 是一個 2 × 1 2\times 1 的截距項向量, A 1 A_1 A 2 A_2 均為 2 × 2 2\times 2 的引數矩陣, ϵ t \epsilon _{ t } 是不與時間相關的擾動項。我假設 ϵ t \epsilon _{ t } 服從分佈 N ( 0 , Σ ) N(0,\Sigma) ,其中 Σ \Sigma 是一個 2 × 2 2\times 2 的協方差矩陣。

設定樣本量為 1000 ,並在 Stata 中生成需要的變數:

 clear all

. set seed 2016

. local T = 1100

. set obs `T'
number of observations (_N) was 0, now 1,100

. gen time = _n

. tsset time
        time variable:  time, 1 to 1100
                delta:  1 unit

. generate y1 = .
(1,100 missing values generated)

. generate y2 = .
(1,100 missing values generated)

. generate eps1 = .
(1,100 missing values generated)

. generate eps2 = .
(1,100 missing values generated)

在第 1-6 行,我設定了隨機種子,將樣本量設為 1000,並且生成了一個時間變數 time 。接下來,我生成了變數 y1y2eps1eps2 來存放觀測序列和擾動項。

2. 設定係數值

我為本文的 VAR(2) 模型選擇瞭如下的引數值:
(2) A 1 = [ 0.6 0.3 0.4 0.2 ] , A 2 = [ 0.2 0.3 0.1 0.1 ] A_{ 1 }= \begin{bmatrix} 0.6 &-0.3 \\ 0.4& 0.2 \end{bmatrix},\quad A_{ 2 }= \begin{bmatrix} 0.2 &0.3 \\ -0.1& 0.1 \end{bmatrix} \tag{2}

. mata:
------------------------------------------------- mata (type end to exit) -----
: mu = (0.1\0.4)

: A1 = (0.6,-0.3\0.4,0.2)

: A2 = (0.2,0.3\-0.1,0.1)

: Sigma = (1,0.5\0.5,1)

: end
-------------------------------------------------------------------------------

在 Mata 中,我分別建立了矩陣 μ \mu A 1 A_1 A 2 A_2 Σ \Sigma 來放置引數值。在模擬得到樣本之前,我檢查了由這些引數值是否會得到一個平穩的 VAR(2) 模型。令:

F = [ A 1 A 2 I 2 0 ] F= \begin{bmatrix} A_{ 1 }&A_{ 2 } \\ I_{ 2 } & 0 \end{bmatrix}

其中 I 2 I_2 2 × 2 2\times 2 的識別矩陣, 0 0 是一個每個元素均為 0 的 2 × 2 2\times 2 矩陣。如果矩陣 F 的所有特徵根均小於 1 ,則此 VAR(2) 過程即為一個平穩過程。下面放上計算特徵根的程式碼:

. mata:
------------------------------------------------- mata (type end to exit) -----
: K = p = 2               // K = number of variables; p = number of lags

: F = J(K*p,K*p,0)

: F[1..2,1..2] = A1

: F[1..2,3..4] = A2

: F[3..4,1..2] = I(K)

: X = L = .

: eigensystem(F,X,L)

: L'
                              1
    +----------------------------+
  1 |                .858715598  |
  2 |  -.217760515 + .32727213i  |
  3 |  -.217760515 - .32727213i  |
  4 |                .376805431  |
    +----------------------------+

: end
------------------------------------------------------------------------------

我根據之前的設定建立了矩陣 F 並使用函式 eigensystem( ) 計算其特徵根。 矩陣 X 中儲存了特徵向量,L 中儲存特徵根。所有 L 中的特徵根均小於 1 ,因此該 VAR(2) 過程是平穩的。在檢驗過是否平穩之後,接下來生成 VAR(2) 模型的擾動項。

3. 由多維正態分佈中生成擾動項

我從分佈 N ( 0 , Σ ) N(0,\Sigma) 中生成兩個隨機正態變數,並將它們分別賦值給變數 eps1eps2

. mata:
------------------------------------------------- mata (type end to exit) -----
: T = strtoreal(st_local("T"))

: u = rnormal(T,2,0,1)*cholesky(Sigma)

: epsmat = .

: st_view(epsmat,.,"eps1 eps2")

: epsmat[1..T,.] = u

: end
-------------------------------------------------------------------------------

我將樣本規模(在 Stata 中定義為一個區域性巨集變數 T )賦值給一個 Mata 數值變數,這步將可以簡化之後的工作。在 Mata 中,我使用了兩個函式:st_local( )strtoreal( ) 來儲存樣本大小。第一個函式可以從 Stata 巨集中獲取文字值,第二個函式則可以將文字值轉變為實數值。

第二行生成了一個 1100 × 2 1100\times 2 的正態擾動項矩陣,其每個元素都服從分佈 N ( 0 , Σ ) N(0,\Sigma) 。我使用函式 st_view( ) 將生成的正態擾動項存放