VAR, SVAR and SVEC Models: Implementation Within R Package vars

The structure of the package vars and its implementation of vector autoregressive, structural vector autoregressive and structural vector error correction models are explained in this paper. In addition to the three cornerstone functions VAR(), SVAR() and SVEC() for estimating such models, functions for diagnostic testing, estimation of a restricted models, prediction, causality analysis, impulse response analysis and forecast error variance decomposition are provided too. It is further possible to convert vector error correction models into their level VAR representation. The different methods and functions are elucidated by employing a macroeconomic data set for Canada. However, the focus in this writing is on the implementation part rather than the usage of the tools at hand.


Introduction
Since the critique of Sims (1980) in the early eighties of the last century, multivariate data analysis in the context of vector autoregressive models (henceforth: VAR) has evolved as a standard instrument in econometrics. 1Because statistical tests are frequently used in determining inter-dependencies and dynamic relationships between variables, this methodology was soon enriched by incorporating non-statistical a priori information.VAR models explain the endogenous variables solely by their own history, apart from deterministic regressors.In contrast, structural vector autoregressive models (henceforth: SVAR) allow the explicit modeling of contemporaneous interdependence between the left-hand side variables.Hence, these types of models try to bypass the shortcomings of VAR models.At the same time as Sims jeopardized the paradigm of multiple structural equation models laid out by the Cowles Foundation in the 1940s and 1950s, Granger (1981) and Engle and Granger (1987) endowed econometricians with a powerful tool for modeling and testing economic relationships, namely, the concept of cointegration.Nowadays these branches of research are unified in the form of vector error correction (henceforth: VECM) and structural vector error correction models (henceforth: SVEC).A thorough theoretical exposition of all these models is provided in the monographs of Lütkepohl (2006), Hendry (1995), Johansen (1995), Hamilton (1994), Banerjee, Dolado, Galbraith, and Hendry (1993).
To the author's knowledge, currently only functions in the base distribution of R and in the CRAN (Comprehensive R Archive Network) packages dse (Gilbert 2000(Gilbert , 1995(Gilbert , 1993) ) and fArma (Würtz 2007) are made available for estimating ARIMA and VARIMA time series models.Although the CRAN package MSBVAR (Brandt and Appleby 2007) provides methods for estimating frequentist and Bayesian vector autoregression (BVAR) models, the methods and functions provided in the package vars try to fill a gap in the econometrics' methods landscape of R by providing the "standard" tools in the context of VAR, SVAR and SVEC analysis.
This article is structured as follows: in the next section the considered models, i.e., VAR, SVAR, VECM and SVEC, are presented.The structure of the package as well as the implemented methods and functions are explained in Section˜3.In the last part, examples of applying the tools contained in vars are exhibited.Finally, a summary and a computational details section conclude this article.

Vector autoregressive models
In its basic form, a VAR consists of a set of K endogenous variables y t = (y 1t , . . ., y kt , . . ., y Kt ) for k = 1, . . .K. The VAR(p)-process is then defined as:2 with A i are (K × K) coefficient matrices for i = 1, . . ., p and u t is a K-dimensional process with E(u t ) = 0 and time invariant positive definite covariance matrix E(u t u t ) = Σ u (white noise).
One important characteristic of a VAR(p)-process is its stability.This means that it generates stationary time series with time invariant means, variances and covariance structure, given sufficient starting values.One can check this by evaluating the characteristic polynomial: If the solution of the above equation has a root for z = 1, then either some or all variables in the VAR(p)-process are integrated of order one, i.e., I(1).It might be the case, that cointegration between the variables does exist.This instance can then be better analyzed in the context of a VECM.
In practice, the stability of an empirical VAR(p)-process can be analyzed by considering the companion form and calculating the eigenvalues of the coefficient matrix.A VAR(p)-process can be written as a VAR(1)-process: with: whereby the dimensions of the stacked vectors ξ t and v t is (KP × 1) and the dimension of the matrix A is (Kp × Kp).If the moduli of the eigenvalues of A are less than one, then the VAR(p)-process is stable.
For a given sample of the endogenous variables y 1 , . . .y T and sufficient presample values y −p+1 , . . ., y 0 , the coefficients of a VAR(p)-process can be estimated efficiently by least-squares applied separately to each of the equations.
Once a VAR(p) model has been estimated, the avenue is wide open for further analysis.A researcher might/should be interested in diagnostic tests, such as testing for the absence of autocorrelation, heteroscedasticity or non-normality in the error process.He might be interested further in causal inference, forecasting and/or diagnosing the empirical model's dynamic behavior, i.e., impulse response functions (henceforth: IRF) and forecast error variance decomposition (henceforth: FEVD).The latter two are based upon the Wold moving average decomposition for stable VAR(p)-processes which is defined as: with Φ 0 = I K and Φ s can be computed recursively according to: whereby A j = 0 for j > p.
Finally, forecasts for horizons h ≥ 1 of an empirical VAR(p)-process can be generated recursively according to: where y T +j|T = y T +j for j ≤ 0. The forecast error covariance matrix is given as: and the matrices Φ i are the empirical coefficient matrices of the Wold moving average representation of a stable VAR(p)-process as shown above.The operator ⊗ is the Kronecker product.

Structural vector autoregressive models
Recall from Section˜2.1 the definition of a VAR(p)-process, in particular Equation˜1.A VAR(p) can be interpreted as a reduced form model.A SVAR model is its structural form and is defined as: It is assumed that the structural errors, ε t , are white noise and the coefficient matrices A * i for i = 1, . . ., p, are structural coefficients that differ in general from their reduced form counterparts.To see this, consider the resulting equation by left-multiplying Equation˜8 with the inverse of A: A SVAR model can be used to identify shocks and trace these out by employing IRA and/or FEVD through imposing restrictions on the matrices A and/or B. Incidentally, though a SVAR model is a structural model, it departs from a reduced form VAR(p) model and only restrictions for A and B can be added.It should be noted that the reduced form residuals can be retrieved from a SVAR model by u t = A −1 Bε t and its variance-covariance matrix by Depending on the imposed restrictions, three types of SVAR models can be distinguished: • B model: A is set to I K (minimum number of restrictions to be imposed for identification is the same as for A model).
• AB model: restrictions can be placed on both matrices (minimum number of restrictions for identification is The parameters are estimated by minimizing the negative of the concentrated log-likelihood function: whereby Σu signifies an estimate of the reduced form variance/covariance matrix for the error process.

Vector error correction models
Reconsider the VAR from Equation˜1: The following vector error correction specifications do exist, which can be estimated with function ca.jo() contained in urca for more details (Pfaff 2006): The Γ i matrices contain the cumulative long-run impacts, hence this VECM specification is signified by "long-run" form.The other specification is given as follows and is commonly employed: with Equation˜14 applies to this specification too.Hence, the Π matrix is the same as in the first specification.However, the Γ i matrices now differ, in the sense that they measure transitory effects.Therefore this specification is signified as "transitory" form.In case of cointegration the matrix Π = αβ is of reduced rank.The dimensions of α and β is K × r and r is the cointegration rank, i.e. how many long-run relationships between the variables y t do exist.The matrix α is the loading matrix and the coefficients of the long-run relationships are contained in β.

Structural vector error correction models
Reconsider the VECM from Equation˜15.It is possible to apply the same reasoning of SVAR models to SVEC models, in particular when the equivalent level-VAR representation of the VECM is used.However, the information contained in the cointegration properties of the variables are thereby not used for identifying restrictions on the structural shocks.Hence, typically a B model is assumed whence a SVEC model is specified and estimated.
whereby u t = Bε t and ε t ∼ N (0, I K ).In order to exploit this information, one considers the Beveridge-Nelson moving average representation of the variables y t if they adhere to the VECM process as in Equation˜15: The variables contained in y t can be decomposed into a part that is integrated of order one and a part that is integrated of order zero.The first term on the right-hand-side of Equation˜18 is referred to the "common trends" of the system and this term drives the system y t .The middle term is integrated of order zero and it is assumed that the infinite sum is bounded, i.e.Ξ * j converge to zero as j → ∞.The initial values are captured by y * 0 .For the modeling of SVEC the interest centers on the common trends in which the long-run effects of shocks are captured.The matrix Ξ is of reduced rank K − r, whereby r is the count of stationary cointegration relationships.The matrix is defined as: Because of its reduced rank only be K − r common trends drive the system.Therefore, by knowing the rank of Π one can then conclude that at most r of the structural errors can have a transitory effect.This implies that at most r columns of Ξ can be set to zero.One can combine the Beveridge-Nelson decomposition with the relationship between the VECM error terms and the structural innovations.The common trends term is then ΞB ∞ t=1 ε t and the long-run effects of the structural innovations are captured by the matrix ΞB.The contemporaneous effects of the structural errors are contained in the matrix B. As in the case of SVAR models of type B one needs for local just-identified SVEC models 1 2 K(K − 1) restrictions.The cointegration structure of the model provides r(K − r) restrictions on the long-run matrix.The remaining restrictions can be placed on either matrix, whereby at least r(r − 1)/2 of them must be imposed directly on the contemporaneous matrix B.

Overview
In Table˜1 the structure of the package vars is outlined.The functions and methods will be addressed briefly here.A more detailed discussion is provided in the following subsections.
When a VAR(p) has been fitted with VAR() one obtains a list object with class attribute varest.
VAR(y, p = 1, type = c("const", "trend", "both", "none"), season = NULL, exogen = NULL, lag.max = NULL, ic = c("AIC", "HQ", "SC", "FPE")) As one can see from Table˜1, basically all relevant investigations can be conducted with the methods and functions made available for this kind of object.Plotting, prediction, forecast error variance decomposition, impulse response analysis, log-likelihood, MA representations and summary are implemented as methods.Diagnostic testing, restricted VAR estimation, stability analysis in the sense of root-stability and empirical fluctuation processes as well as causality analysis are implemented as functions.Some of the methods have their own print and plot methods.Furthermore, extractor methods for obtaining the residuals, the fitted values and the estimated coefficients do exist.
In Section˜2.2 it has been argued that a VAR can be viewed as particular SVAR model.The function SVAR() requires therefore an object of class varest.The default estimation method of a SVAR model is a scoring algorithm, as proposed by Amisano and Giannini (1997).
The restrictions for the A, B or A and B matrices have to be provided in explicit form as arguments Amat and/or Bmat.Alternatively, the different SVAR types can be estimated by directly minimizing the negative log-likelihood.The latter is used if the estimation method has then to be set to estmethod = "direct".
Structural vector error correction models can be estimated with the function SVEC().
Finally, objects of formal class ca.jo generated with function ca.jo() contained in the package urca can be converted to their level VAR representation with the function vec2var().The resultant object has class attribute vec2var and the analytical tools are likewise applicable as in the case of objects with class attribute varest.

Cornerstone functions
The function for estimating a VAR(p) model is VAR().It consists of seven arguments.A data matrix or an object that can be coerced to it has to be provided for y.The lag-order has to be submitted as integer for p.In general, this lag-order is unknown and VAR() offers the possibility to automatically determine an appropriate lag inclusion.This is achieved by by setting lag.max to an upper bound integer value and ic to a desired information criteria.Possible options are Akaike (ic = "AIC"), Hannan-Quinn (ic = "HQ"), Schwarz (ic = "SC") or the forecast prediction error (ic = "FPE").The calculations are based upon the same sample size.That is lag.maxvalues are used for each of the estimated models as starting values. 3The type of deterministic regressors to include into the VAR(p) is set by the argument type.Possible values are to include a constant, a trend, both or none deterministic regressors.In addition, the inclusion of centered seasonal dummy variables can be achieved by setting season to the seasonal frequency of the data (e.g., for quarterly data: season = 4).Finally, exogenous variables, like intervention dummies, can be included by providing a matrix object for exogen.
The summary method returns -aside of descriptive information about the estimated VARthe estimated equations as well as the variance/covariance and the correlation matrix of the residuals.It is further possible to report summary results for selected equations only.This is achieved by using the function's argument equation which expects a character vector with the names of the desired endogenous variables.The plot method displays for each equation in a VAR a diagram of fit, a residual plot, the auto-correlation and partial auto-correlation function of the residuals in a single layout.If the plot method is called interactively, the user is requested to enter <RETURN> for commencing to the next plot.It is further possible to plot the results for a subset of endogenous variables only.This is achieved by using the argument name in the plot method.The appearance of the plot can be adjusted to ones liking by setting the relevant arguments of the plot method.
A SVAR model is estimated with the function SVAR().An object with class attribute varest has to be provided as argument for x.The structural parameters are estimated either by a scoring algorithm (the default) or by direct minimization of the negative log-likelihood function.Whether an A, B or AB˜model will be estimated, is dependent on the setting for Amat and Bmat.If a restriction matrix for Amat with dimension (K × K) is provided and the argument Bmat is left NULL, an A model will be estimated.In this case Bmat is set internally to an identity matrix I K .Alternatively, if only a matrix object for Bmat is provided and Amat is left unchanged, then a B model will be estimated and internally Amat is set to an identity matrix I K .Finally, if matrix objects for both arguments are provided, then an AB model will be estimated.In all cases, the matrix elements to be estimated are marked by NA entries at the relevant positions.The user has the option to provide starting values for the unknown coefficients by providing a vector object for the argument start.If start is left NULL, then 0.1 will be used as the starting value for all coefficients.The arguments max.iter, conv.crit and maxls can be used for tuning the scoring algorithm.The maximum number of iterations is controlled by max.iter, the convergence criterion is set by conv.crit and the maximum step length is set by maxls.A likelihood ratio test is computed per default for overidentified systems.This default setting can be offset by lrtest = FALSE.If a just-identified has been set-up, a warning is issued that an over-identification test cannot be computed in case of lrtest = TRUE.The ellipsis argument (...) is passed to optim() in case of direct optimization.
The returned object of function SVAR() is a list with class attribute svarest.Dependent on the chosen model and if the argument hessian = TRUE has been set in case of estmethod = "direct", the list elements A, Ase, B, Bse contain the estimated coefficient matrices with the numerical standard errors.The element LRIM does contain the longrun impact matrix in case a SVAR of type Blanchard & Quah is estimated with function BQ(), otherwise this element is NULL (see Blanchard and Quah 1989).The list element Sigma.U is the variance-covariance matrix of the reduced form residuals times 100, i.e., Σ U = A −1 BB A −1 × 100.The list element LR is an object with class attribute htest, holding the likelihood ratio over-identification test.The element opt is the returned object from function optim() in case estmethod = "direct" has been used.The remaining five list items are the vector of starting values, the SVAR model type, the varest object, the number of iterations and the call to SVAR().
A SVEC model is estimated with the function SVEC().The supplied object for the argument x must be of formal class ca.jo.The restrictions on the long-run and short-run structural coefficient matrices must be provided as arguments LR and SR, respectively.These matrices have either zero or NA entries as their elements.It is further necessary to specify the cointegration rank of the estimated VECM via the argument r.Likewise to SVAR(), the arguments start, max.iter, conv.crit and maxls can be used for tuning the scoring algorithm.The argument lrtest applies likewise as in SVAR().Finally, the logical flag boot can be used for calculating standard errors by applying the bootstrap method to the SVEC.The count of repetition is set by the argument runs.The returned list object from SVEC() and its associated methods will be bespoken in Section 4, where a SVEC model is specified for the Canadian labor market.
Finally, with function vec2var() a VECM (i.e., an object of formal class ca.jo, generated by the function ca.jo() contained in the package urca) is transformed into its level-VAR representation.Aside of this argument the function requires the cointegrating rank as this information is needed for the transformation (see Equations˜12-16).The print method does return the coefficient values, first for the lagged endogenous variables, next for the deterministic regressors.

Diagnostic testing
In the package vars functions for diagnostic testing are arch.test(),normality.test(),serial.test()and stability().The former three functions return a list object with class attribute varcheck for which plot and print method exist.The plots -one for each equation -include a residual plot, an empirical distribution plot and the ACF and PACF of the residuals and their squares.The plot method offers additional arguments for adjusting its appearance.
The implemented tests for heteroscedasticity are the univariate and multivariate ARCH test (see Engle 1982;Hamilton 1994;Lütkepohl 2006).The multivariate ARCH-LM test is based on the following regression (the univariate test can be considered as special case of the exhibition below and is skipped): whereby v t assigns a spherical error process and vech is the column-stacking operator for symmetric matrices that stacks the columns from the main diagonal on downward.The dimension of β 0 is 1 2 K(K + 1) and for the coefficient matrices B i with i = 1, . . ., q, 1 2 K(K + 1) × 1 2 K(K + 1).The null hypothesis is: H 0 := B 1 = B 2 = . . .= B q = 0 and the alternative is: H 1 : B 1 = 0 ∩ B 2 = 0 ∩ . . .∩ B q = 0.The test statistic is defined as: with and Ω assigns the covariance matrix of the above defined regression model.This test statistic is distributed as χ 2 (qK 2 (K + 1) 2 /4).
The The returned tests have class attribute htest, hence the print method for these kind of objects is implicitly used as print method for objects of class varcheck.This applies likewise to the functions normality.test()and serial.test(),which will be bespoken next.
The Jarque-Bera normality tests for univariate and multivariate series are implemented and applied to the residuals of a VAR(p) as well as separate tests for multivariate skewness and kurtosis (see Bera andJarque 1980, 1981;Jarque and Bera 1987;Lütkepohl 2006).The univariate versions of the Jarque-Bera test are applied to the residuals of each equation.A multivariate version of this test can be computed by using the residuals that are standardized by a Choleski decomposition of the variance-covariance matrix for the centered residuals.Please note, that in this case the test result is dependent upon the ordering of the variables.The test statistics for the multivariate case are defined as: whereby s 2 3 and s 2 4 are computed according to: with b 1 and b 2 are the third and fourth non-central moment vectors of the standardized residuals ûs t = P − (û t − ūt ) and P is a lower triangular matrix with positive diagonal such that P P = Σu , i.e., the Choleski decomposition of the residual covariance matrix.The test statistic JB mv is distributed as χ 2 (2K) and the multivariate skewness, s 2 3 , and kurtosis test, s 2 4 are distributed as χ 2 (K).

normality.test(x, multivariate.only = TRUE)
The matrix of residuals is the first element in the returned list object.Edgerton and Shukur (1999).
The Portmanteau statistic is defined as: with Ĉi = 1 T Σ T t=i+1 ût û t−i .The test statistic has an approximate χ 2 (K 2 h − n * ) distribution, and n * is the number of coefficients excluding deterministic terms of a VAR(p).The limiting distribution is only valid for h tending to infinity at a suitable rate with growing sample size.Hence, the trade-off is between a decent approximation to the χ 2 distribution and a loss in power of the test, when h is chosen too large.The small sample adjustment of the test statistic is given as: and is computed if type = "PT.adjusted" is set.
The Breusch-Godfrey LM statistic is based upon the following auxiliary regressions (see Breusch 1978;Godfrey 1978): The null hypothesis is: 0 and correspondingly the alternative hypothesis is of the form H 1 : ∃B i = 0 f or i = 1, 2, . . ., h.The test statistic is defined as: where ΣR and Σe assign the residual covariance matrix of the restricted and unrestricted model, respectively.The test statistic LM h is distributed as χ 2 (hK 2 ).Edgerton and Shukur (1999) proposed a small sample correction, which is defined as: ), whereby n is the number of regressors in the original system and m = Kh.The modified test statistic is distributed as F (hK 2 , int(N r − q)).This test statistic is returned if type = "ES" has been used.serial.test(x,lags.pt= 16, lags.bg= 5, type = c("PT.asymptotic","PT.adjusted", "BG", "ES")) The test statistics are returned in the list element serial and have class attribute htest.
Per default the asymptotic Portmanteau test is returned.The residuals are contained in the first list element.
The function stability() returns a list object with class attribute varstabil.The function itself is just a wrapper for the function efp() contained in the package strucchange (see Zeileis, Leisch, Hornik, and Kleiber 2002, for a detailed exposition of the package's capabilities).The first element of the returned list object is itself a list of objects with class attribute efp.
Hence, the plot method for objects of class varstabil just calls the plot method for objects of class efp.

Prediction, impulse responses and forecast error decomposition
A predict method for objects with class attribute varest or vec2var is available.The n.ahead forecasts are computed recursively for the estimated VAR(p)-process and a value for the forecast confidence interval can be provided too.Its default value is 0.95.The confidence interval is inferred from the empirical forecast error covariance matrix.
predict(object, ..., n.ahead = 10, ci = 0.95, dumvar = NULL) The predict method returns a list with class attribute varprd.The forecasts are contained as a list in its first element signified as fcst.The second entry are the endogenous variables themselves.The last element is the submitted model object to predict.The print method returns the forecasts with upper and lower confidence levels, if applicable.The plot method draws the time series plots, whereby the start of the out-of-sample period is marked by a dashed vertical line.If this method is called interactively, the user is requested to browse through the graphs for each variable by hitting the <RETURN> key.For visualizing forecasts (i.e., objects with class attribute varprd) fan charts can be generated by the function fanchart() (see Britton, Fisher, and Whitley 1998).If the functional argument color is not set, the colors are taken from the gray color scheme.Likewise, if no confidence levels are supplied, then the fan charts are produced for confidence levels ranging from 0.1 to 0.9 with step size 0.1.
fanchart(x, colors = NULL, cis = NULL, names = NULL, main = NULL, ylab = NULL, xlab = NULL, col.y = NULL, nc, plot.type= c("multiple", "single"), mar = par("mar"), oma = par("oma"), ...) The impulse response analysis is based upon the Wold moving average representation of a VAR(p)-process (see Equations˜5 and 6 above).It is used to investigate the dynamic interactions between the endogenous variables.The (i, j)th coefficients of the matrices Φ s are thereby interpreted as the expected response of variable y i,t+s to a unit change in variable y jt .These effects can be accumulated through time, s = 1, 2, . .., and hence one would obtain the simulated impact of a unit change in variable j to the variable i at time s.Aside of these impulse response coefficients, it is often conceivable to use orthogonal impulse responses as an alternative.This is the case, if the underlying shocks are less likely to occur in isolation, but when contemporaneous correlations between the components of the error process u t exist, i.e., the off-diagonal elements of Σ u are non-zero.The orthogonal impulse responses are derived from a Choleski decomposition of the error variance-covariance matrix: Σ u = P P with P being a lower triangular.The moving average representation can then be transformed to: with ε t = P −1 u t and Ψ i = Φ i P for i = 0, 1, 2, . . .and Ψ 0 = P .Incidentally, because the matrix P is lower triangular, it follows that only a shock in the first variable of a VAR(p)process does exert an influence on all the remaining ones and that the second and following variables cannot have a direct impact on y 1t .Please note, that a different ordering of the variables might produce different outcomes with respect to the impulse responses.The nonuniqueness of the impulse responses can be circumvented by analyzing a set of endogenous variables in the SVAR framework.
Impulse response analysis has been implemented as a method for objects with class attribute of either varest, svarest, svecest or vec2var.These methods are utilizing the methods Phi and Psi, where applicable.
irf(x, impulse = NULL, response = NULL, n.ahead = 10, ortho = TRUE, cumulative = FALSE, boot = TRUE, ci = 0.95, runs = 100, seed = NULL, ...) The impulse variables are set as a character vector impulse and the responses are provided likewise in the argument response.If either one is unset, then all variables are considered as impulses or responses, respectively.The default length of the impulse responses is set to 10 via argument n.ahead.The computation of orthogonal and/or cumulative impulse responses is controlled by the logical switches ortho and cumulative, respectively.Finally, confidence bands can be returned by setting boot = TRUE (default).The preset values are to run 100 replications and return 95% confidence bands.It is at the user's leisure to specify a seed for the random number generator.The standard percentile interval is calculated as , where s * γ/2 and s * (1−γ)/2 are the γ/2 and (1 − γ)/2 quantiles of the estimated bootstrapped impulse response coefficients Φ * or Ψ * (see Efron and Tibshirani 1993).The irf method returns a list object with class attribute varirf for which print and plot methods do exist.Likewise to the plot method of objects with class attribute varprd, the user is requested to browse through the graphs for each variable by hitting the <RETURN> key, whence the method is called interactively.The appearance of the plots can be adjusted.
The forecast error variance decomposition is based upon the orthogonal impulse response coefficient matrices Ψ n .The FEVD allows the user to analyze the contribution of variable j to the h-step forecast error variance of variable k.If the element-wise squared orthogonal impulse responses are divided by the variance of the forecast error variance, σ 2 k (h), the resultant is a percentage figure.The fevd method is available for conducting FEVD.Methods for objects of classes varest, svarest, svecest and vec2var do exist.Aside of the object itself, the argument n.ahead can be specified; its default value is 10.-summary(ur.df(diff(Canada[, "prod"]), type = "drift", + lags = 1)) R> adf2 According to the AIC and FPE the optimal lag number is p = 3, whereas the HQ criterion indicates p = 2 and the SC criterion indicates an optimal lag length of p = 1.They estimated for all three lag orders a VAR including a constant and a trend as deterministic regressors and conducted diagnostic tests with respect to the residuals.In the R˜code example below, the relevant commands are exhibited for the VAR(1) model.First, the variables have to be reordered in the same sequence as in Breitung et˜al. (2004).This step is necessary, because otherwise the results of the multivariate Jarque-Bera test, in which a Choleski decomposition is employed, would differ slightly from the reported ones in Breitung et˜al. (2004).In the R˜code lines below the estimation of the VAR(1) as well as the summary output and the diagram of fit for equation "e" is shown.
R> Canada <-Canada[, c("prod", "e", "U", "rw")] R> p1ct <-VAR(Canada, p = 1, type = "both") R> p1ct   Given the diagnostic test results the authors concluded that a VAR(1)-specification might be too restrictive.They argued further, that although some of the stability tests do indicate deviations from parameter constancy, the time-invariant specification of the VAR(2) and VAR(3) model will be maintained as tentative candidates for the following cointegration analysis.
The authors estimated a VECM whereby a deterministic trend has been included in the cointegration relation.The estimation of these models as well as the statistical inference with respect to the cointegration rank can be swiftly accomplished with the function ca.jo().
Although the following R˜code examples are using functions contained in the package urca, it is however beneficial to reproduce these results for two reasons: the interplay between the functions contained in package urca and vars is exhibited and it provides an understanding of the then following SVEC specification.
R> summary(ca.jo(Canada,type = "trace", ecdet = "trend", K = 3, + spec = "transitory")) R> summary(ca.jo(Canada,type = "trace", ecdet = "trend", K = 2, + spec = "transitory")) The outcome of the trace tests is provided in Table˜4.These results do indicate one cointegration relationship.The reported critical values differ slightly from the ones that are reported in Table 4.3 of Breitung et˜al. (2004).The authors used the values that are contained in Johansen (1995), whereas the values from Osterwald-Lenum (1992) are used in the function ca.jo().In the R˜code snippet below the VECM is re-estimated with this restriction and a normalization of the long-run relationship with respect to real wages.The results are shown in Table˜5.
R> vecm <-ca.jo(Canada[,c("rw", "prod", "e", "U")], type = "trace", + ecdet = "trend", K = 3, spec = "transitory") R> vecm.r1<-cajorls(vecm, r = 1) For a just identified SVEC model of type B one needs 1 2 K(K − 1) = 6 linear independent restrictions.It is further reasoned from the Beveridge-Nelson decomposition that there are k * = r(K − r) = 3 shocks with permanent effects and only one shock that exerts a temporary effect, due to r = 1.Because the cointegration relation is interpreted as a stationary wagesetting relation, the temporary shock is associated with the wage shock variable.Hence, the four entries in the last column of the long-run impact matrix ΞB are set to zero.Because this matrix is of reduced rank, only k * r = 3 linear independent restrictions are imposed thereby.
It is therefore necessary to set 1 2 k * (k * − 1) = 3 additional elements to zero.The authors assumed constant returns of scale and therefore productivity is only driven by output shocks.This reasoning implies zero coefficients in the first row of the long-run matrix for the variables employment, unemployment and real wages, hence the elements ΞB 1,j for j = 2, 3, 4 are set to zero.Because ΞB 1,4 has already been set to zero, only two additional restrictions have been added.The last restriction is imposed on the element B 4,2 .Here, it is assumed that labor demand shocks do not exert an immediate effect on real wages.
In the R˜code example below the matrix objects LR and SR are set up accordingly and the just-identified SVEC is estimated with function SVEC().In the call to the function SVEC() the argument boot = TRUE has been employed such that bootstrapped standard errors and hence t˜statistics can be computed for the structural long-run and contemporaneous coefficients.
The authors investigated further if labor supply shocks do have no long-run impact on unemployment.This hypothesis is mirrored by setting ΞB 3,3 = 0.Because one more zero restriction has been added to the long-run impact matrix, the SVEC model is now over-identified.The validity of this over-identification restriction can be tested with a LR test.In the R˜code example below first the additional restriction has been set and then the SVEC is re-estimated by using the upate method.The result of the LR test is contained in the returned list as named element LRover.The value of the test statistic is 6.07 and the p˜value of this χ 2 (1)-distributed variable is 0.014.Therefore, the null hypothesis that shocks to the labor supply do not exert a long-run effect on unemployment has to be rejected for a significance level of 5%.
In order to investigate the dynamic effects on unemployment, the authors applied an impulse response analysis.The impulse response analysis shows the effects of the different shocks, i.e., output, labor demand, labor supply and wage, to unemployment.In the R˜code example below the irf method for objects with class attribute svecest is employed and the argument boot = TRUE has been set such that confidence bands around the impulse response trajectories can be calculated.R> svec.irf<-irf(svec, response = "U", n.ahead = 48, boot = TRUE) R> plot(svec.irf) The outcome of the IRA is exhibited in Figure˜5.
In a final step, a forecast error variance decomposition is conducted with respect to unemployment.This is achieved by applying the fevd method to the object with class attribute svecest.
R> fevd.U <-fevd(svec, n.ahead = 48)$U The authors report only the values for selected quarters.These results are displayed in Table˜8.

Summary
In this paper it has been described how the different functions and methods contained in the package vars are designed and offer the researcher a fairly easy to use environment for conducting VAR, SVAR and SVEC analysis.This is primarily achieved through methods for impulse response functions, forecast error variance decomposition, prediction as well as tools for diagnostic testing, determination of a suitable lag length for the model and stability / causality analysis.
The package vars complements the package urca in the sense that most of the above described tools are available for VECM that can be swiftly transformed to their level-VAR representation whence the cointegrating rank has been determined.

Computational details
The package's code is purely written in R (R Development Core Team 2008) and S3-classes with methods have been utilized.It is shipped with a NAMESPACE and a ChangeLog file.It has dependencies to MASS (Venables and Ripley 2002), strucchange (Zeileis et˜al. 2002) and urca (Pfaff 2006).R itself as well as these packages can be obtained from CRAN at http://CRAN.R-project.org/.Furthermore, daily builds of package vars are provided on R-Forge (see http://R-Forge.R-project.org/projects/vars/).It has been published under GPL version˜2 or newer.The results used in this paper were obtained using R 2.7.0 with packages vars 1.3-8, strucchange 1.3-3, urca 1.1-6 and MASS 7.2-42.

Figure 1 :
Figure 1: Canadian labor market time series.

Figure 5 :
Figure 5: Responses of unemployment to economic shocks with a 95% bootstrap confidence interval.

Table 1 :
Structure of package vars.
default is to compute the multivariate test only.If multivariate.only= FALSE, the univariate tests are computed too.In this case, the returned list object from arch.test() has three elements.The first element is the matrix of residuals.The second, signified by arch.uni, is a list object itself and holds the univariate test results for each of the series.The multivariate test result is contained in the third list element, signified by arch.mul.
arch.test(x, lags.single= 16, lags.multi= 5, multivariate.only= TRUE) The function's default is to compute the multivariate test statistics only.The univariate tests are returned if multivariate.only= FALSE is set.Similar to the returned list object of function arch.test() for this case, the univariate versions of the Jarque-Bera test are applied to the residuals of each equation and are contained in the second list element signified by jb.uni.The multivariate version of the test as well as multivariate tests for skewness and kurtosis are contained in the third list element signified by jb.mul.For testing the lack of serial correlation in the residuals of a VAR(p), a Portmanteau test and the Breusch-Godfrey LM test are implemented in the function serial.test().For both tests small sample modifications can be calculated too, whereby the modification for the LM test has been introduced by

Table 3 :
The resulting graphic is displayed in Figure˜2.Next, it is shown how the diagnostic tests are conducted for the VAR(1) model.The results of all diagnostic tests, i.e., for the VAR(1), VAR(2) and VAR(3) model, are provided in Table˜3 and the graphics of the varcheck object arch1 for the employment equation and the OLS-CUSUM tests for the VAR(1) model are shown in Figures˜3 and 4, respectively.Diagnostic tests of VAR(p) specifications for Canadian data.

Table 4 :
Johansen cointegration tests for Canadian system.

Table 5 :
Cointegration vector and loading parameters (with t˜statistics in parentheses).

Table 6 :
Estimated coefficients of the contemporaneous impact matrix (with t˜statistics in parentheses).

Table 7 :
Estimated coefficients of the long-run impact matrix (with t˜statistics in parentheses).

Table 8 :
Forecast error variance decomposition of Canadian unemployment.