MULTIVARIATE STOCHASTIC PROCESSES

DEFINITION

A multivariate time series (x₁, ... , x_t, ... , x_n) with x_t = [ x_1t, ... , x_kt] is a finite realization of a multivariate stochastic process {X_t}_t∈[0,∞] where X_t is a random vector defined on an appropriate probability space (Ω, Β_R, P), with expected value equal to:
E[x_t]=μ_t=[ μ_1t, ... , μ_kt]
and autocovariance matrix equal to:
E[x_t, x_t-j]= E[(x_t-μ_t) (x_t-j-μ_t-j)'] ----> for j ≠ 0
E[x_t, x_t-j]= Var(x_t-μ_t) ----> for j = 0

PROPERTIES

Strong Stationarity:

A stochastic process is said to be strongly stationary if the distribution function of the process remains invariant to transformations.
Weak Stationarity:

A stochastic process is said to be weakly stationary if the expected value and autocovariance matrix do not vary as t varies.
Ergodicity:

^-1

Invertibility:

_t-1

_t-2

MULTIVARIATE WHITE NOISE

A white-noise process is a weakly stationary process that has zero mean and is uncorrelated over time. A multivariate white noise is defined as the stochastic process {u_t = [u_1t,...,u_kt]}_t∈[0,∞] consisting of a sequence of k-dimensional random variables that are uncorrelated when referring to different time instants.
In other words, the fact that a multivariate process is white noise excludes correlation between each element of the process and the past history of the entire process, but does not exclude that there may be correlation between contemporary elements.

VARMA MODELS

A class of stationary stochastic processes, applicable to multiple phenomena considered jointly, is constituted by VARMA(p,q) models. These models represent a generalization of ARMA models used in the univariate context.
A VARMA(p,q) model can be expressed in the form:

X_t = c + Φ₁X_t-1 + ... + Φ_pX_t-p + u_t + Θ₁u_t-1 + ... + Θ_qu_t-q

u_t is the multivariate white noise process;
c indicates a kx1 vector of deterministic constants;
Φ_j is a kxk matrix of autoregressive coefficients;
Θ_i are the coefficients for the moving average part.

VARMA models are rarely used in empirical applications for two reasons. First, the completeness of such a model creates significant problems regarding the precision of estimating the process parameters. Moreover, since in many circumstances the autoregressive part exhaustively explains the time series under examination (the moving average part has negligible relevance in the phenomenon's evolution, as it is due to transitory shocks that vanish over time), it is often preferable to neglect the moving average component, thus obtaining a particular case of the VARMA model, called VAR (vector autoregressive models).

VAR MODELS

DEFINITION

We now consider a particular case of the VARMA(p,q) model where all coefficients of the moving average part Θ_j are set equal to zero. The resulting model is called the VAR(p) model (endowed with the invertibility property) whose analytical expression is:

(X_t - μ) = Φ₁(X_t-1 - μ) + ... + Φ_p(X_t-p - μ) + u_t

equivalently:

(I_k - Φ₁L- ... - Φ_pL^p)(X_t - μ) = u_t

or in compact form:

Φ(L)(X_t - μ) = u_t

where ut is a white noise. If we set p = 1 we obtain the VAR(1) model:

(X_t - μ) = Φ1(X_{t- 1} - μ) + u_t

VAR IN COMPANION FORM

The companion form representation allows us to represent any VAR(p) model as a VAR(1). Starting from the following model:

(X_t - μ) = Φ1(X_{t- 1} - μ) + u_t

We construct a system of equations where the first is the one just written and the others are:

(X_t-1 - μ) = (X_t-1 - μ)
...
...
...
(X_t-p+1 - μ) = (X_t-p+1 - μ)

we obtain the companion form as:

ε_t=Fε_t-1+v_t

STATIONARITY OF THE VAR

Given a process VAR(1) it can be defined as weakly stationary if the autovalues λ_i of the matrix Φ₁ are inside the unit circle. The same reasoning can be extended to a VAR(p) process written in companion form as:
ε_t=Fε_t-1+v_t

where:

ε_t is a (kp x 1) vector => [X_t - μ, X_t-1 - μ, ... , X_t-p+1 - μ]
F is a matrix (kp x kp) containing in the first row the vector of the Φ vertically concatenated with the Kronecker product of an identity matrix of order k by an identity matrix of order p-1 concatenated with a zero vector of width p-1
v_t vt is a (kp x 1) vector, whose first component is a vector of k white noises.

For a process X_t to be defined as stationary, all eigenvalues of the matrix F must be inside the unit circle.

FROM VAR(∞) AL VAR(1)

If the eigenvalues λ_i of the matrix Φ₁ are inside the unit circle (have modulus less than 1) then:

Φ₁^j -_{_j->∞}-> 0,
∑^∞_j=0 Φ₁^j = (I_k - Φ₁)^-1 exists and is finite;

furthermore, the fact that the eigenvalues of the matrix Φ₁ are inside the unit circle guarantees that Φ₁^j -_{_j->∞}-> 0 and that ∑_j=0^∞ Φ₁^j = (I_k - Φ₁)^-1, which implies the absolute summability of this same series. With these two conditions, the VAR(1) process can be written as VMA(∞), and this automatically guarantees stationarity.

ESTIMATION OF A VAR MODEL

The procedure for identifying the best VAR model is not immediate, but occurs through the comparison of different models.
Model estimation is carried out through a maximum likelihood approach, and since the variables are assumed to be normally distributed, maximum likelihood estimates and Ordinary Least Squares coincide.
The first step is to establish the order p of the model. This choice can be made through an appropriate balance between number of parameters and significance of the likelihood reduction.
In practice this is done using loss functions, minimising a function that is increasing with respect to the number of parameters and decreasing with respect to the likelihood.
let us consider a VAR(p) model:

(X_t - μ) = Φ₁(X_t-1 - μ) + ... + Φ_p(X_t-p - μ) + u_t

dove u_t ∼ WN(0 , Σ)

With the ML method we want to obtain estimates for Φ and Σ; however, since the variables are correlated, we cannot write the likelihood function as a product of the marginals, so it is more convenient to consider the conditional likelihood function:

L^C(Φ₁, ... , Φ_p ; Σ) = ƒ_{_{_{X_t+n, ... ,X_t+1}|_{X_t=x_t, ... , X_t-p+1=x_t-p+1}}}(X_t+n, ... , X_t+1)=

ⁿ

_i=1

_{_{_{X_t+i}|_{x_t, ... , x_t-p+1}}}

_t+i

It remains to identify the ƒ_{_{_{X_t+i}|_{x_t, ... , x_t-p+1}}}(X_t+i), analizzando the process X_t+i (which is a VAR(p)) and knowing that u_t ∼ WN ∼ N(0 , Σ) we can arrive at the conclusion that:

X_t+i|x_t, ... , x_t-p+1 ∼ N(Π^T Z_t+i , Σ)

where: X_t+i = c + Φ₁X_t+i + ... + Φ_p+iX_t-p+1
and therefore:
Π^T = [c , Φ₁ , ... + Φ_p+i]
Z_t+i = [1 , X_t+i , ... , X_t-p+1]^T

Knowing that the conditional distribution of X_t+i|(given its past) ∼ N(Π^T Z_t+i , Σ) we can proceed with ML estimation:

(Π^ML)^{^T} = [∑ⁿ_i=1 (X_t+i Z^T_t+i)][∑ⁿ_i=1 (X_t+i Z^T_t+i)]^-1
Σ^ML = n^-1 ∑ⁿ_i=1 {[X_t+i - (Π^ML)^{^T}Z_t+i][X_t+i - (Π^ML)^{^T}Z_t+i]^T}

RECURSIVE LAG EXCLUSION TEST

The recursive lag exclusion test is based on the likelihood ratio and aims to identify the best lag order p for the VAR(p) model:
H₀: p = p₀
H₁: p = p₁ (t.c. p₁ > p₀)

We consider:

λ(x) = ^{L^C(Π₀^ML , Σ₀^ML)} / _{L^C(Π₁^ML , Σ₁^ML)}

According to Wald's theorem we know that:

W(x) = -2log(λ(x)) = -2[l^C(Π₀^ML , Σ₀^ML) - l^C(Π₁^ML , Σ₁^ML)]

where: W(x) ∼ χ²_{_p₁-p₀}

This means that the rejection region will be:
R={x tc W(x)>q_{(α, χ²_{_p₁-p₀} )} }
Where q_{(α, χ²_{_p₁-p₀}} ) is the quantile di ordine α di a χ²_{_p₁-p₀}.

Another way to identify the lag order p is to adopt methods that minimise functions that are increasing in the parameters and decreasing with respect to the likelihood; these functions are called loss functions. The main ones used in this context are the Akaike, Hannan-Quinn, and Schwarz functions:

AIC(p) = -2 n^-1 l(θ , p) + 2 n^-1 p k² -> higher p

HQ(p) = -2 n^-1 l(θ , p) + 2 n^-1 log(n^-1 log(n)) p k²

AIC(p) = -2 n^-1 l(θ , p) + 2 n^-1 log(n) p k² -> lower p

A third method, although not very practical, can be to evaluate multiple models based on the correlations with the residuals, keeping in mind that a good model should have approximately 95% of the residuals uncorrelated.

#### VAR Models
    data(Canada) # macroeconomic data on Canada
    summary(Canada)

##        and              prod             rw              U         
##  Min.   :928.6   Min.   :401.3   Min.   :386.1   Min.   : 6.700  
##  1st Qu.:935.4   1st Qu.:404.8   1st Qu.:423.9   1st Qu.: 7.782  
##  Median :946.0   Median :406.5   Median :444.4   Median : 9.450  
##  Mean   :944.3   Mean   :407.8   Mean   :440.8   Mean   : 9.321  
##  3rd Qu.:950.0   3rd Qu.:410.7   3rd Qu.:461.1   3rd Qu.:10.607  
##  Max.   :961.8   Max.   :418.0   Max.   :470.0   Max.   :12.770

    # Lag order selection
    # Lag order selection with 4 criteria 
    VARselect(Canada, lag.max = 12, type="const")

## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      3      1      1      3 
## 
## $criteria
##                   1            2            3            4            5
## AIC(n) -6.655129808 -6.814854295 -6.868744361 -6.666586167 -6.460700426
## HQ(n)  -6.403366607 -6.361680533 -6.214160038 -5.810591283 -5.403294982
## SC(n)  -6.022722553 -5.676521235 -5.224485498 -4.516401499 -3.804589954
## FPE(n)  0.001288555  0.001103138  0.001056564  0.001319817  0.001676682
##                   6            7            8            9          10
## AIC(n) -6.363384957 -6.292604661 -6.279055217 -6.201794055 -6.17085150
## HQ(n)  -5.104568952 -4.832378095 -4.617418090 -4.338746367 -4.10639325
## SC(n)  -3.201348681 -2.624642580 -2.105167332 -1.521980366 -0.98511201
## FPE(n)  0.001944191  0.002244786  0.002518411  0.003123173  0.00387657
##                  11           12
## AIC(n) -6.229490861 -6.447761913
## HQ(n)  -3.963622052 -3.980482542
## SC(n)  -0.537825564 -0.250170811
## FPE(n)  0.004681472  0.005242163

     # Estimation of the VAR model with 3 lags
    m1<-VAR(Canada, p = 3, type = "const")

    summary(m1)

## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: e, prod, rw, U 
## Deterministic variables: const 
## Sample size: 81 
## Log Likelihood: -150.609 
## Roots of the characteristic polynomial:
## 1.004 0.9283 0.9283 0.7437 0.7437 0.6043 0.6043 0.5355 0.5355 0.2258 0.2258 0.1607
## Call:
## VAR(y = Canada, p = 3, type = "const")
## 
## 
## Estimation results for equation e: 
## ================================== 
## and = e.l1 + prod.l1 + rw.l1 + U.l1 + e.l2 + prod.l2 + rw.l2 + U.l2 + e.l3 + prod.l3 + rw.l3 + U.l3 + const 
## 
##           Estimate Std. Error t value Pr(>|t|)    
## e.l1       1.75274    0.15082  11.622  < 2e-16 ***
## prod.l1    0.16962    0.06228   2.723 0.008204 ** 
## rw.l1     -0.08260    0.05277  -1.565 0.122180    
## U.l1       0.09952    0.19747   0.504 0.615915    
## e.l2      -1.18385    0.23517  -5.034 3.75e-06 ***
## prod.l2   -0.10574    0.09425  -1.122 0.265858    
## rw.l2     -0.02439    0.06957  -0.351 0.727032    
## U.l2      -0.05077    0.24534  -0.207 0.836667    
## e.l3       0.58725    0.16431   3.574 0.000652 ***
## prod.l3    0.01054    0.06384   0.165 0.869371    
## rw.l3      0.03824    0.05365   0.713 0.478450    
## U.l3       0.34139    0.20530   1.663 0.100938    
## const   -150.68737   61.00889  -2.470 0.016029 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.3399 on 68 degrees of freedom
## Multiple R-Squared: 0.9988,	Adjusted R-squared: 0.9985 
## F-statistic:  4554 on 12 and 68 DF,  p-value: < 2.2e-16 
## 
## 
## Estimation results for equation prod: 
## ===================================== 
## prod = e.l1 + prod.l1 + rw.l1 + U.l1 + e.l2 + prod.l2 + rw.l2 + U.l2 + e.l3 + prod.l3 + rw.l3 + U.l3 + const 
## 
##           Estimate Std. Error t value Pr(>|t|)    
## e.l1      -0.14880    0.28913  -0.515   0.6085    
## prod.l1    1.14799    0.11940   9.615 2.65e-14 ***
## rw.l1      0.02359    0.10117   0.233   0.8163    
## U.l1      -0.65814    0.37857  -1.739   0.0866 .  
## e.l2      -0.18165    0.45083  -0.403   0.6883    
## prod.l2   -0.19627    0.18069  -1.086   0.2812    
## rw.l2     -0.20337    0.13337  -1.525   0.1319    
## U.l2       0.82237    0.47034   1.748   0.0849 .  
## e.l3       0.57495    0.31499   1.825   0.0723 .  
## prod.l3    0.04415    0.12239   0.361   0.7194    
## rw.l3      0.09337    0.10285   0.908   0.3672    
## U.l3       0.40078    0.39357   1.018   0.3121    
## const   -195.86985  116.95813  -1.675   0.0986 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.6515 on 68 degrees of freedom
## Multiple R-Squared:  0.98,	Adjusted R-squared: 0.9765 
## F-statistic: 277.5 on 12 and 68 DF,  p-value: < 2.2e-16 
## 
## 
## Estimation results for equation rw: 
## =================================== 
## rw = e.l1 + prod.l1 + rw.l1 + U.l1 + e.l2 + prod.l2 + rw.l2 + U.l2 + e.l3 + prod.l3 + rw.l3 + U.l3 + const 
## 
##           Estimate Std. Error t value Pr(>|t|)    
## e.l1    -4.716e-01  3.373e-01  -1.398    0.167    
## prod.l1 -6.500e-02  1.393e-01  -0.467    0.642    
## rw.l1    9.091e-01  1.180e-01   7.702 7.63e-11 ***
## U.l1    -7.941e-04  4.417e-01  -0.002    0.999    
## e.l2     6.667e-01  5.260e-01   1.268    0.209    
## prod.l2 -2.164e-01  2.108e-01  -1.027    0.308    
## rw.l2   -1.457e-01  1.556e-01  -0.936    0.353    
## U.l2    -3.014e-01  5.487e-01  -0.549    0.585    
## e.l3    -1.289e-01  3.675e-01  -0.351    0.727    
## prod.l3  2.140e-01  1.428e-01   1.498    0.139    
## rw.l3    1.902e-01  1.200e-01   1.585    0.118    
## U.l3     1.506e-01  4.592e-01   0.328    0.744    
## const   -1.167e+01  1.365e+02  -0.086    0.932    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.7601 on 68 degrees of freedom
## Multiple R-Squared: 0.9989,	Adjusted R-squared: 0.9987 
## F-statistic:  5239 on 12 and 68 DF,  p-value: < 2.2e-16 
## 
## 
## Estimation results for equation U: 
## ================================== 
## U = e.l1 + prod.l1 + rw.l1 + U.l1 + e.l2 + prod.l2 + rw.l2 + U.l2 + e.l3 + prod.l3 + rw.l3 + U.l3 + const 
## 
##          Estimate Std. Error t value Pr(>|t|)    
## e.l1     -0.61773    0.12508  -4.939 5.39e-06 ***
## prod.l1  -0.09778    0.05165  -1.893 0.062614 .  
## rw.l1     0.01455    0.04377   0.332 0.740601    
## U.l1      0.65976    0.16378   4.028 0.000144 ***
## e.l2      0.51811    0.19504   2.656 0.009830 ** 
## prod.l2   0.08799    0.07817   1.126 0.264279    
## rw.l2     0.06993    0.05770   1.212 0.229700    
## U.l2     -0.08099    0.20348  -0.398 0.691865    
## e.l3     -0.03006    0.13627  -0.221 0.826069    
## prod.l3  -0.01092    0.05295  -0.206 0.837180    
## rw.l3    -0.03909    0.04450  -0.879 0.382733    
## U.l3      0.06684    0.17027   0.393 0.695858    
## const   114.36732   50.59802   2.260 0.027008 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.2819 on 68 degrees of freedom
## Multiple R-Squared: 0.9736,	Adjusted R-squared: 0.969 
## F-statistic: 209.2 on 12 and 68 DF,  p-value: < 2.2e-16 
## 
## 
## 
## Covariance matrix of residuals:
##             and     prod       rw        U
## and     0.11550 -0.03161 -0.03681 -0.07034
## prod -0.03161  0.42449  0.05589  0.01494
## rw   -0.03681  0.05589  0.57780  0.03660
## U    -0.07034  0.01494  0.03660  0.07945
## 
## Correlation matrix of residuals:
##            and     prod      rw        U
## and     1.0000 -0.14276 -0.1425 -0.73426
## prod -0.1428  1.00000  0.1129  0.08136
## rw   -0.1425  0.11286  1.0000  0.17084
## U    -0.7343  0.08136  0.1708  1.00000

    roots(m1) # roots of the VAR model

##  [1] 1.0038607 0.9282722 0.9282722 0.7436701 0.7436701 0.6043160 0.6043160
##  [8] 0.5354599 0.5354599 0.2257507 0.2257507 0.1606576

    roots(m1,modulus = F) # roots of the VAR model as complex numbers

##  [1]  1.0038607+0.0000000i  0.9264902+0.0574900i  0.9264902-0.0574900i
##  [4]  0.7048399+0.2371623i  0.7048399-0.2371623i  0.1092774+0.5943536i
##  [7]  0.1092774-0.5943536i -0.1181711+0.5222575i -0.1181711-0.5222575i
## [10]  0.1907349+0.1207625i  0.1907349-0.1207625i -0.1606576+0.0000000i

    residuals(m1)

##               and        prod           rw           U
## 1   0.396121511  0.11840048  0.659939985 -0.43738364
## 2   0.249065788  0.71742949  0.431450548  0.15319901
## 3   0.177175520  0.06623809  1.384352350 -0.09961548
## 4  -0.096351638 -0.93100815 -1.237251416  0.05339949
## 5  -0.034436534  0.06846126  2.075215983  0.12308452
## 6  -0.447350640  0.49219100  0.115836902 -0.14146982
## 7  -0.609808396 -0.40496047  0.293517252  0.53160459
## 8  -0.615737798  1.16668022 -0.643492743  0.61839378
## 9   0.006032917  0.39733305  0.120832350 -0.04799867
## 10  0.167250285  0.65813337 -1.438020699 -0.40110084
## 11 -0.067311190 -0.72145313  0.394209375  0.14627210
## 12  0.039902231 -0.74312679 -0.284892821  0.16797709
## 13 -0.575903232  0.31992989  1.147899743  0.16852248
## 14 -0.177616108  0.48710456  0.128979113  0.17964754
## 15  0.036238860  0.33431182 -0.290204411  0.13443535
## 16  0.019971532 -0.60775077 -0.168174576  0.09141991
## 17 -0.385827418  0.79944851  0.423511106  0.25676295
## 18 -0.204161553 -0.08778962  0.171456986  0.08972733
## 19  0.316398779 -0.59706989 -0.431137686 -0.24631910
## 20 -0.180180306 -0.28901892 -0.448742807 -0.04331506
## 21  0.120755711  0.21324333  0.682417761  0.19992276
## 22  0.195147342 -1.82297173 -0.423433462 -0.15443427
## 23  0.170930944 -0.09592562 -0.961001661 -0.01989833
## 24 -0.169875846 -0.39776715 -1.861015184  0.01160187
## 25  0.177713435 -0.59720908 -0.003871509 -0.18328396
## 26  0.138432688  1.25522727 -0.372485011  0.00607115
## 27  0.522406319 -0.47422733 -1.090729372 -0.29221673
## 28 -0.023033217  0.48798380 -0.509418297 -0.04023790
## 29  0.607301887  0.35781522 -0.063729533 -0.36376680
## 30 -0.114757003  0.41345841  0.145481150  0.08017616
## 31  0.042105446  0.27619096 -0.701423817 -0.16808739
## 32 -0.030203576  0.50730605 -0.892996698 -0.07211574
## 33  0.373238302  0.50677007 -0.096810513 -0.25955166
## 34  0.322595594 -0.54798254  0.845465989 -0.03033181
## 35 -0.526183816  0.54261376 -0.437697818  0.13801490
## 36  0.610300275 -0.33401206 -0.533955938 -0.65741436
## 37 -0.006201963 -0.72940223  0.377727739 -0.05148297
## 38  0.393125175 -0.35703835 -0.277260277 -0.18815652
## 39  0.146653356 -0.06573476 -0.121862691 -0.22434903
## 40 -0.062053825 -0.41348199 -1.511185348  0.15788602
## 41 -0.580460268 -0.37317205  0.483649371  0.48488056
## 42 -0.292238455 -1.98327299  1.148691052  0.39484445
## 43  1.074654073  0.29325269  0.106699570 -0.56772324
## 44 -0.192627916 -0.53614979  0.044547943  0.01772929
## 45 -0.094668664  0.28829404 -0.019961688 -0.18211791
## 46 -0.305803217 -0.17647904 -0.141635494 -0.04669920
## 47  0.073986217 -0.20275269  0.066039910 -0.03145155
## 48 -0.084687188  0.17669476 -0.028737430  0.21583533
## 49  0.093636224 -0.04885420  0.187185399  0.31932022
## 50  0.342289135 -0.48998088  0.797749434 -0.49631186
## 51 -0.304197228  0.39808535  0.074324443  0.72742843
## 52  0.214058552 -0.40404990 -0.220003411 -0.10973062
## 53 -0.100527096 -0.02646415 -0.124656561  0.10080570
## 54 -0.201497545  1.14485721  0.745536655 -0.16039754
## 55  0.276082101 -0.04585363 -0.010666951 -0.13412455
## 56  0.006836110  0.14169369  0.676564181 -0.17719865
## 57  0.019940849 -0.20230280  0.693085332 -0.06006858
## 58  0.092266176 -0.30079762 -0.203785433  0.01363578
## 59 -0.479901663  0.21355976 -0.127031468 -0.03626645
## 60  0.251008227  0.04232691  0.118932383 -0.26098548
## 61 -0.123598672  0.02038972 -0.138870691 -0.34236644
## 62 -0.036743625 -0.75955610 -0.138363811  0.09763768
## 63  0.031496469  0.39870925  1.022046660 -0.18355803
## 64 -0.030422843  1.14512515  1.233092101  0.23061948
## 65 -0.120871818  0.87571369  0.311475352  0.14241668
## 66  0.666117698  0.09228093 -0.450304496 -0.52819034
## 67 -0.237500396 -0.42224627  0.067040762  0.17470641
## 68  0.363035033  0.18137274 -1.100665633 -0.15400239
## 69 -0.362961805  0.26810366  1.271041966  0.33129146
## 70 -0.091926400 -0.29819929  0.473261543  0.08544725
## 71  0.352029192 -0.73971295  0.483579939 -0.25626999
## 72  0.194373788 -0.82628942 -1.080886819 -0.06163559
## 73  0.178071734  0.73743633 -0.548341524  0.02461780
## 74 -0.122256032  0.56425442 -0.471141542  0.25750303
## 75 -0.008935260  0.02182941  0.015663746  0.29060378
## 76 -0.159509669  0.66632598  0.328149401  0.03111804
## 77 -0.079783374 -0.06483946 -0.333436999 -0.22821044
## 78 -0.122557116  0.48111926  1.097083694  0.11731999
## 79 -0.662671358  0.30624382 -0.278514408  0.31481300
## 80 -0.324176024  0.39010860  0.360424935  0.34626384
## 81 -0.011227787 -0.91514419 -0.986361458  0.11288173

    # representation of the roots on the unit circle
    plot(c(-1, 1), c(-1, 1), type = "n")
    radius <- 1
    theta <- seq(0, 2 * pi, length = 200)
    lines(x = radius * cos(theta), y = radius * sin(theta))
    abline(h=0); abline(v=0)
    points(roots(m1,modulus = F),col=4,pch=16)

USES OF THE VAR

Forecasting
The best predictor is the one that minimises the mean squared error (MSE), i.e. the expected value of the random variable X_t+1 conditioned on the information available at time t. Therefore the forecast at time t+1 is given by:
estimate -> x_{_t+1|_{ℑ_t}} ≡ E[X_t+1|ℑ_t]

Impulse Response
The impulse response function is useful for recording the reactions of the system to sudden shocks with respect to each of the k variables considered in the process. The impulse response function also serves to evaluate the persistence characteristics of a given process.

Forecast Variance Decomposition
The error made in predicting the evolution of the VAR, h periods ahead, is given by the expression:
MSE[E[X_(t+h)|ℑ_t]] = Σ + Θ₁ Σ Θ'₁ + ... + Θ_h-1 Σ Θ'_h-1
dove Σ = E[u_t u'_t]

Granger Causality Analysis
The correlation between two variables implies the existence of a link, but we do not know the type of cause-effect relationship that connects them.
Granger studied this relationship; in particular, if X_t and Y_t are two time series, he states that Y_t does not Granger-cause X_t if the forecast of X_t+s based on information from X at time t is not worse than the same forecast obtained based on information from both x_t and Y_t

#### Uses of the VAR

    #Forecasting
    library(vars)
    data(Canada)
    summary(Canada)

##        and              prod             rw              U         
##  Min.   :928.6   Min.   :401.3   Min.   :386.1   Min.   : 6.700  
##  1st Qu.:935.4   1st Qu.:404.8   1st Qu.:423.9   1st Qu.: 7.782  
##  Median :946.0   Median :406.5   Median :444.4   Median : 9.450  
##  Mean   :944.3   Mean   :407.8   Mean   :440.8   Mean   : 9.321  
##  3rd Qu.:950.0   3rd Qu.:410.7   3rd Qu.:461.1   3rd Qu.:10.607  
##  Max.   :961.8   Max.   :418.0   Max.   :470.0   Max.   :12.770

    m1<-VAR(Canada, p = 3, type = "const")
    predict(m1, n.ahead = 5, ci = 0.95)

## $e
##          fcst    lower    upper        CI
## [1,] 962.8355 962.1694 963.5016 0.6661063
## [2,] 964.1396 962.8205 965.4587 1.3190726
## [3,] 965.5223 963.6664 967.3782 1.8559305
## [4,] 966.8774 964.6137 969.1411 2.2636929
## [5,] 968.1325 965.5151 970.7499 2.6173821
## 
## $prod
##          fcst    lower    upper       CI
## [1,] 417.4415 416.1645 418.7185 1.276970
## [2,] 418.0955 416.1377 420.0533 1.957804
## [3,] 418.8073 416.3501 421.2645 2.457200
## [4,] 419.3108 416.4527 422.1689 2.858117
## [5,] 419.6458 416.4004 422.8911 3.245365
## 
## $rw
##          fcst    lower    upper       CI
## [1,] 470.0419 468.5521 471.5318 1.489829
## [2,] 470.8222 468.7613 472.8832 2.060919
## [3,] 471.2779 468.9004 473.6554 2.377508
## [4,] 471.7219 469.0327 474.4110 2.689163
## [5,] 472.3777 469.3794 475.3759 2.998213
## 
## $U
##          fcst    lower    upper        CI
## [1,] 6.484233 5.931795 7.036672 0.5524385
## [2,] 5.831105 4.922149 6.740061 0.9089557
## [3,] 5.166632 3.918600 6.414664 1.2480320
## [4,] 4.576521 3.065484 6.087559 1.5110372
## [5,] 4.051549 2.350973 5.752125 1.7005762

    # Impulse Response (Dynamic Analysis)
    imp<-irf(m1, impulse = c("e","prod", "rw", "U"),
             response = c("e","prod", "rw", "U"),
             boot = FALSE,n.ahead = 40)
    plot(imp)

    plot(imp$irf$U[,1],type="l",ylim=c(-0.4,1))
    abline(h=0)
    lines(imp$irf$U[,2],lty=2,col="blue")
    lines(imp$irf$U[,3],lty=3,col="green")
    lines(imp$irf$U[,4],lty=4,col="gray")

    #Forecast Variance Decomposition
    fevd(m1, n.ahead = 5)

## $e
##              and       prod          rw            U
## [1,] 1.0000000 0.00000000 0.000000000 0.0000000000
## [2,] 0.9679295 0.02335231 0.007926938 0.0007912516
## [3,] 0.8901016 0.07017553 0.039079270 0.0006436471
## [4,] 0.7767949 0.13515474 0.083448990 0.0046014017
## [5,] 0.6413324 0.20621412 0.128869199 0.0235842576
## 
## $prod
##                and      prod           rw          U
## [1,] 0.020380831 0.9796192 0.000000e+00 0.00000000
## [2,] 0.009233159 0.9750352 2.297056e-05 0.01570867
## [3,] 0.010597113 0.9666723 9.394910e-03 0.01333570
## [4,] 0.019496778 0.9415248 2.689223e-02 0.01208623
## [5,] 0.037266067 0.9025632 3.977396e-02 0.02039678
## 
## $rw
##               and        prod        rw            U
## [1,] 0.02029883 0.008738054 0.9709631 0.000000e+00
## [2,] 0.06827676 0.005031703 0.9266915 2.063692e-08
## [3,] 0.07114374 0.037957529 0.8885843 2.314384e-03
## [4,] 0.05593555 0.072967669 0.8689615 2.135299e-03
## [5,] 0.05387566 0.088640298 0.8542818 3.202191e-03
## 
## $U
##              and        prod          rw          U
## [1,] 0.5391325 0.000562155 0.004824606 0.45548072
## [2,] 0.7333695 0.020725917 0.004418956 0.24148560
## [3,] 0.7815156 0.045194065 0.033812485 0.13947783
## [4,] 0.7327095 0.084688479 0.084300593 0.09830144
## [5,] 0.6426027 0.143944586 0.135094921 0.07835783

    # Example of Granger Causality Analysis
    data(ChickEgg)
    summary(ChickEgg)

##     chicken            egg      
##  Min.   :364584   Min.   :3081  
##  1st Qu.:387658   1st Qu.:5008  
##  Median :403819   Median :5380  
##  Mean   :419504   Mean   :4986  
##  3rd Qu.:433773   3rd Qu.:5530  
##  Max.   :582197   Max.   :5836

    dataGr<-as.data.frame(ChickEgg)
    par(mfrow=c(2,1))
    plot.ts(dataGr$chicken)
    plot.ts(dataGr$egg)

    dchick <- diff(dataGr$chicken)
    degg <- diff(dataGr$egg)
    plot.ts(dchick)
    plot.ts(degg)

     # p-value below 0.05 therefore there is a difference between the two models
    grangertest(dchick ~ degg, order=4)

## Granger causality test
## 
## Model 1: dchick ~ Lags(dchick, 1:4) + Lags(degg, 1:4)
## Model 2: dchick ~ Lags(dchick, 1:4)
##   Res.Df Df      F   Pr(>F)   
## 1     40                      
## 2     44 -4 4.1762 0.006414 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

     # p-value above 0.05 no difference between the two models
    grangertest(degg ~ dchick, order=4)

## Granger causality test
## 
## Model 1: degg ~ Lags(degg, 1:4) + Lags(dchick, 1:4)
## Model 2: degg ~ Lags(degg, 1:4)
##   Res.Df Df      F Pr(>F)
## 1     40                 
## 2     44 -4 0.2817 0.8881

USES OF THE VAR

Once the estimates of the VAR(p) coefficients are obtained, we may find ourselves in a non-stationarity situation if some eigenvalue λ_i of the matrix F (matrix of the VAR(p) coefficients in companion form) is found outside or on the unit circle.
The situation is more uncertain if some eigenvalue is found very close to 1; in this case we should test stationarity with a non-stationarity test such as the Dickey-Fuller test, since in case of non-stationarity all results regarding estimator consistency would fail; in case of non-stationarity, regression can be used.

A problem that arises if we decide to switch to regression but obtain coefficients equal to 0 (which would imply stationarity) is that of residual correlation.
Two cases must be distinguished:

a correlation is called spurious (or nonsense correlation) if it links two phenomena that have no causal connection, that have nothing in common;

Y US Defence Spending, 1971-1990

X Population in South Africa, 1971-1990

Y = -368.99 + 0.0179 X
(-11.34) (16.75)

R^2 = 0.940, D/W = 0.4069, F = 280.69

Corr = .9694

it dice invece indiretta when due variables X and Y are correlate perchè in reality correlate entrambe a a variable Z.

A classic example is the one linking the variable X=number of victims in a fire to the variable Y=number of firefighters employed to extinguish it. A positive correlation is easily found: more firefighters, more victims; it would therefore be "logical", but completely absurd, to think that employing fewer firefighters would result in fewer victims. Obviously X and Y depend on the third variable Z=size of the fire.

As a solution to these problems, we can differentiate to remove the stochastic trend when two variables are correlated.

Python in Practice

Below we simulate and estimate a VAR(1) model using Python with fake data.

1. Simulating a VAR(1) Process

import numpy as np
from statsmodels.tsa.api import VAR

np.random.seed(42)
n = 300

# VAR(1) coefficient matrix
# y1_t = 0.5*y1_{t-1} + 0.1*y2_{t-1} + e1
# y2_t = 0.2*y1_{t-1} + 0.6*y2_{t-1} + e2
A = np.array([[0.5, 0.1],
              [0.2, 0.6]])

data = np.zeros((n, 2))
for t in range(1, n):
    data[t] = A @ data[t-1] + np.random.normal(0, 1, 2)

print("Sample means: Y1={}, Y2={}".format(round(data[:,0].mean(), 3), round(data[:,1].mean(), 3)))
print("Sample stds:  Y1={}, Y2={}".format(round(data[:,0].std(), 3), round(data[:,1].std(), 3)))

# Output:
# Sample means: Y1=0.156, Y2=0.247
# Sample stds:  Y1=1.491, Y2=1.563

2. Model Selection and Estimation

model = VAR(data)
results = model.fit(maxlags=5, ic='aic')
print("Selected lag order (AIC):", results.k_ar)
print("\nEstimated coefficients (lag 1):")
print(results.coefs[0].round(3))
print("\nTrue coefficients:")
print(A)

# Output:
# Selected lag order (AIC): 1
#
# Estimated coefficients (lag 1):
# [[0.487 0.101]
#  [0.185 0.613]]
#
# True coefficients:
# [[0.5 0.1]
#  [0.2 0.6]]

3. Granger Causality Test

granger = results.test_causality('y1', 'y2', kind='f')
print("Granger causality: Y2 -> Y1")
print(f"  F-statistic: {granger.test_statistic:.4f}")
print(f"  p-value: {granger.pvalue:.4f}")
print(f"  Causal: {'Yes' if granger.pvalue < 0.05 else 'No'}")

# Output:
# Granger causality: Y2 -> Y1
#   F-statistic: 3.0214
#   p-value: 0.0832
#   Causal: No (marginal - true coeff is only 0.1)

4. Impulse Response Function and Forecast

irf = results.irf(20)
print("IRF of Y1 to a shock in Y1 (first 5 periods):")
print(irf.irfs[:5, 0, 0].round(3))

forecast = results.forecast(data[-results.k_ar:], steps=10)
print("\n10-step forecast:")
print("  Y1:", forecast[:, 0].round(3))
print("  Y2:", forecast[:, 1].round(3))

# Output:
# IRF of Y1 to a shock in Y1 (first 5 periods):
# [1.    0.487 0.256 0.141 0.08 ]
#
# 10-step forecast:
#   Y1: [ 0.674  0.38   0.222  0.13   0.076  0.044  0.026  0.015  0.009  0.005]
#   Y2: [-0.261 -0.009  0.068  0.055  0.039  0.025  0.016  0.01   0.006  0.004]