（一）R语言vars包来进行VAR、SVAR和SVECM的建模

  VAR模型

   VAR（y, p=1, type=c("const","trend", "both", "none"), season=NULL， exogen=NULL, lag.max=NULL,ic=c("AIC", "HQ", "SC", "FPE")) y是一个数据矩阵；p为滞后阶数；type为回归的类型；season为季节数据频数的设置（如果是季节数据，则season=4）；exogen为外生变量矩阵；lag.max为最大滞后阶数；ic为信息准则。

  SVAR模型

    SVAR(x,estmethod=c("scoring", "direct"), Amat=NULL, Bmat=NULL, start=NULL,max.iter=100, conv.crit=1e-07, maxls=1, lrtest=TRUE,...) x为数据矩阵；estmethod为由得分算法（默认）或负对数似然函数直接最小化估计；Amat和Bmat为A和B的结构系数矩阵；start为未知系数的初始值；conv.crit为收敛值；maxls为最大步长；lrtest为逻辑值，用来确定计算过度识别系统的似然比率。

  SVEC模型

   SVEC(x, LR=NULL, SR=NULL,r=1, start=NULL, max.iter=100, conv.crit=1e-07, maxls=1,lrtest=TRUE, boot=FALSE, runs=100) x必须为ca.jo类型格式的，也就是说必须由设定协整方程后再构建SVEC；LR和SR分别为设定长期和短期的结构系数矩阵；r为协整方程个数。

  诊断检验

    vars包中进行诊断检验的函数为arch.test(), normality.test(),serial.test()和stability()

  脉冲响应和预测方差误差分解

  脉冲响应函数

   irf(x, impulse=NULL,response=NULL, n.ahead=10, ortho=TRUE, cumulative=FALSE, boot=TRUE,ci=0.95, runs=100,seed=NULL, ...)

  方差分解

  fevd(x,n.ahead=10,...)

 （二）案例分析

 加拿大的一个宏观经济数据

prod为GDP，e为就业的对数，U为失业率，rw为实际工资水平的对数。  library(vars) data(Canada) summary(Canada)       e           prod           rw            U         Min.  :928.6   Min.  :401.3   Min.  :386.1   Min.  : 6.700    1stQu.:935.4   1st Qu.:404.8   1stQu.: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    3rdQu.:950.0   3rd Qu.:410.7   3rdQu.:461.1   3rd Qu.:10.607   Max.  :961.8   Max.  :418.0   Max.  :470.0   Max.  :12.770   plot(Canada,nc=2,xlab="") adf1<-summary(ur.df(Canada[,"prod"],type="trend",lag=2)) adf1 ###############################################  # Augmented Dickey-FullerTest Unit Root Test #  ###############################################  Test regressiontrend  Call: lm(formula = z.diff ~z.lag.1 + 1 + tt + z.diff.lag) Residuals:     Min      1Q   Median      3Q    Max  -2.19924 -0.38994 0.04294  0.41914 1.71660  Coefficients:            EstimateStd. Error t value Pr(>|t|)   (Intercept) 30.415228 15.309403   1.987  0.0506 . z.lag.1    -0.075791   0.038134 -1.988   0.0505. tt         0.013896   0.006422  2.164   0.0336* z.diff.lag1 0.284866   0.114359  2.491   0.0149* z.diff.lag2 0.080019   0.116090  0.689   0.4927  --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’1 Residual standard error:0.6851 on 76 degrees of freedom Multiple R-squared: 0.1354, Adjusted R-squared: 0.08993  F-statistic: 2.976 on 4 and76 DF,  p-value: 0.02438 Value of test-statistic is:-1.9875 2.3 2.3817  Critical values for teststatistics:       1pct  5pct10pct tau3 -4.04 -3.45-3.15 phi2  6.50 4.88  4.16 phi3  8.73 6.49  5.47 adf2<-summary(ur.df(diff(Canada[,"prod"]),type="drift",lag=1)) adf2 ###############################################  # Augmented Dickey-FullerTest Unit Root Test #  ###############################################  Test regressiondrift  Call: lm(formula = z.diff ~z.lag.1 + 1 + z.diff.lag) Residuals:     Min      1Q   Median      3Q    Max  -2.05124 -0.39530 0.07819  0.41109 1.75129  Coefficients:           Estimate Std. Error t valuePr(>|t|)     (Intercept) 0.11534   0.08029   1.437   0.155    z.lag.1    -0.68893   0.13350  -5.160 1.83e-06*** z.diff.lag -0.04274   0.11275  -0.379   0.706    --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’1 Residual standard error:0.6971 on 78 degrees of freedom Multiple R-squared: 0.3615, Adjusted R-squared: 0.3451  F-statistic: 22.08 on 2 and78 DF,  p-value: 2.526e-08 Value of test-statistic is:-5.1604 13.3184  Critical values for teststatistics:       1pct  5pct10pct tau2 -3.51 -2.89-2.58 phi1  6.70 4.71  3.86 VARselect(Canada,lag.max=8,type="both") $selection AIC(n) HQ(n)  SC(n)FPE(n)      3     2     1    3  $criteria                1          2          3          4          5          6          7 AIC(n) -6.272579064-6.636669705 -6.771176872 -6.634609210 -6.398132246 -6.307704843-6.070727259 HQ(n) -5.978429449 -6.146420347 -6.084827770-5.752160366 -5.319583658 -5.033056512-4.599979185 SC(n) -5.536558009 -5.409967947 -5.053794411-4.426546046 -3.699388378 -3.118280272-2.390621985 FPE(n) 0.001889842  0.001319462 0.001166019  0.001363175 0.001782055  0.002044202 0.002768551               8 AIC(n)-6.06159685 HQ(n) -4.39474903 SC(n) -1.89081087 FPE(n) 0.00306012   从以上结果可知，AIC和FPE选择3阶之后，SC选择1阶。 Canada<-Canada[,c("prod","e","U","rw")] p1ct<-VAR(Canada,p=1,type="both") p1ct VAR EstimationResults: =======================  Estimated coefficients forequation prod:  =========================================  Call: prod = prod.l1 + e.l1 + U.l1+ rw.l1 + const + trend     prod.l1       e.l1       U.l1      rw.l1      const     trend   0.96313671 0.01291155  0.21108918-0.03909399 16.24340747 0.04613085  Estimated coefficients forequation e:  ======================================  Call: e = prod.l1 + e.l1 + U.l1 +rw.l1 + const + trend       prod.l1        e.l1         U.l1       rw.l1        const       trend    0.19465028   1.23892283   0.62301475   -0.06776277-278.76121138  -0.04066045  Estimated coefficients forequation U:  ======================================  Call: U = prod.l1 + e.l1 + U.l1 +rw.l1 + const + trend      prod.l1        e.l1       U.l1       rw.l1       const      trend   -0.12319201 -0.24844234   0.39158002  0.06580819 259.98200967  0.03451663  Estimated coefficients forequation rw:  =======================================  Call: rw = prod.l1 + e.l1 + U.l1 +rw.l1 + const + trend      prod.l1        e.l1       U.l1       rw.l1       const      trend   -0.22308744 -0.05104397  -0.36863956  0.94890946 163.02453066  0.07142229  summary(p1ct,equation="e")#选择就业方程 VAR EstimationResults: =========================  Endogenous variables: prod,e, U, rw  Deterministic variables:both  Sample size:83  Log Likelihood:-207.525  Roots of the characteristicpolynomial: 0.9504 0.9504 0.90450.7513 Call: VAR(y = Canada, p = 1, type= "both") Estimation results forequation e:  ==================================  e = prod.l1 + e.l1 + U.l1 +rw.l1 + const + trend          Estimate Std. Error t value Pr(>|t|)    prod.l1   0.19465   0.03612   5.389 7.49e-07*** e.l1      1.23892   0.08632  14.353 < 2e-16 *** U.l1      0.62301   0.16927   3.681 0.000430*** rw.l1    -0.06776   0.02828  -2.396 0.018991 *  const  -278.76121   75.18295  -3.7080.000392 *** trend    -0.04066   0.01970  -2.064 0.042378 *  --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’1 Residual standard error:0.4701 on 77 degrees of freedom Multiple R-Squared: 0.9975,Adjusted R-squared: 0.9973  F-statistic: 6088 on 5 and 77 DF,  p-value:< 2.2e-16  Covariance matrix ofresiduals:         prod       e       U      rw prod 0.469517  0.06767 -0.04128 0.002141 e    0.067667  0.22096 -0.13200-0.082793 U   -0.041280 -0.13200  0.12161 0.063738 rw   0.002141 -0.08279  0.06374 0.593174 Correlation matrix ofresiduals:         prod      e      U       rw prod 1.000000  0.2101 -0.1728 0.004057 e    0.210085  1.0000 -0.8052-0.228688 U   -0.172753 -0.8052  1.0000 0.237307 rw   0.004057 -0.2287  0.2373 1.000000 plot(p1ct,names="e")#画出就业对数的拟合、残差图和相关图  ser11<-serial.test(p1ct,lags.pt=16,type="PT.asymptotic") ser11$serial Portmanteau Test(asymptotic) data: Residuals of VAR object p1ct Chi-squared = 233.5, df =240, p-value = 0.606 norm1<-normality.test(p1ct) norm1$jb.mul $JB JB-Test(multivariate) data: Residuals of VAR object p1ct Chi-squared = 9.9189, df =8, p-value = 0.2708 $Skewness Skewness only(multivariate) data: Residuals of VAR object p1ct Chi-squared = 6.356, df = 4,p-value = 0.1741 $Kurtosis Kurtosis only(multivariate) data: Residuals of VAR object p1ct Chi-squared = 3.5629, df =4, p-value = 0.4684 arch1<-arch.test(p1ct,lags.multi=5) arch1$arch.mul ARCH(multivariate) data: Residuals of VAR object p1ct Chi-squared = 570.14, df =500, p-value = 0.01606 plot(arch1,names="e")  plot(stability(p1ct),nc=2)  

summary(ca.jo(Canada,type="trace",ecdet="trend",K=3,spec="transitory")) #协整检验 ######################  # Johansen-Procedure#  ######################  Test type: trace statistic ,with linear trend incointegration  Eigenvalues(lambda): [1] 4.505013e-011.962777e-01 1.676668e-01 4.647108e-022.632104e-17 Values of teststatistic andcritical values of test:         test 10pct  5pct 1pct r <= 3 | 3.85 10.49 12.25 16.26 r <= 2 | 18.72 22.7625.32 30.45 r <= 1 | 36.42 39.0642.44 48.45 r = 0  |84.92 59.14 62.99 70.05 Eigenvectors, normalised tofirst column: (These are the cointegrationrelations)            prod.l1      e.l1      U.l1     rw.l1   trend.l1 prod.l1  1.00000000  1.0000000  1.0000000 1.0000000  1.0000000 e.l1    -0.02385143  1.2946681-2.8831559  4.2418087 -8.2903941 U.l1     3.16874549 3.4036732 -7.4261514  6.8413561-12.5578436 rw.l1    1.83528156 -0.3330945 1.3978789 -0.1393999  2.4466500 trend.l1 -1.30156097-0.2302803 -0.5093218 -1.5925918  0.2831079 Weights W: (This is the loadingmatrix)           prod.l1       e.l1       U.l1       rw.l1     trend.l1 prod.d -0.006535281-0.02763446 -0.070975296 -0.014754352 1.077469e-11 e.d   -0.008503348  0.11414011-0.008156659  0.003988051 7.400296e-12 U.d   -0.004718574 -0.06154306 0.020719431 -0.006557248-4.663893e-12 rw.d  -0.046213350 -0.14579644 -0.016945105  0.011896044 6.952035e-12 summary(ca.jo(Canada,type="trace",ecdet="trend",K=2,spec="transitory")) ######################  # Johansen-Procedure#  ######################  Test type: trace statistic ,with linear trend incointegration  Eigenvalues(lambda): [1] 4.483918e-012.323995e-01 1.313250e-01 4.877895e-029.508809e-17 Values of teststatistic andcritical values of test:         test 10pct  5pct 1pct r <= 3 | 4.10 10.49 12.25 16.26 r <= 2 | 15.65 22.7625.32 30.45 r <= 1 | 37.33 39.0642.44 48.45 r = 0  |86.12 59.14 62.99 70.05 Eigenvectors, normalised tofirst column: (These are the cointegrationrelations)           prod.l1       e.l1      U.l1     rw.l1  trend.l1 prod.l1  1.0000000   1.0000000  1.00000000 1.000000  1.000000 e.l1     2.7132129 -6.3190324  0.49616472 16.333916-10.368563 U.l1     8.8369211 -15.2682881 1.48062661 25.774259 -16.048489 rw.l1   -0.3716323   3.1817254-0.04085215 -2.546391   4.927457 trend.l1 -0.4177976 -0.9335588 -0.26592659 -3.413555 -1.753060 Weights W: (This is the loadingmatrix)           prod.l1       e.l1      U.l1       rw.l1     trend.l1 prod.d 0.023155644 -0.02832697 -0.10914770 -0.006295988-4.784701e-13 e.d    0.005602438 -0.01739149 0.08679396 -0.001019323-4.385546e-13 U.d   -0.019277135  0.01381763-0.03696147 -0.002276871 4.919886e-13 rw.d  -0.084618968 -0.02739056 -0.07798404  0.003985020-1.032315e-13 vecm<-ca.jo(Canada[,c("rw","prod","e","U")],type="trace",ecdet="trend",K=3,spec="transitory") vecm.r1<-cajorls(vecm,r=1)#VEC模型 vecm.r1 $rlm Call: lm(formula =substitute(form1), data = data.mat) Coefficients:         rw.d      prod.d     e.d      U.d      ect1     -0.084815 -0.011994  -0.015606 -0.008660 constant 55.469125   8.274808 10.331308  5.687832 rw.dl1   -0.012082   0.004707 -0.078491  0.017263 prod.dl1 -0.074493   0.234441  0.200953 -0.138916 e.dl1    -0.634084  -0.246544  0.821558 -0.646846 U.dl1     0.063137 -0.979868   0.003379 -0.191125 rw.dl2   -0.157388  -0.190264 -0.095835  0.080354 prod.dl2 -0.251940  -0.029520  0.048273 -0.002909 e.dl2     0.081197 -0.580473  -0.459693 -0.019741 U.dl2    -0.230009  -0.128101 -0.103415 -0.262685 $beta              ect1 rw.l1    1.00000000 prod.l1  0.54487553 e.l1    -0.01299605 U.l1     1.72657188 trend.l1-0.70918872 vecm<-ca.jo(Canada[,c("prod","e","U","rw")],type="trace",ecdet="trend",K=3,spec="transitory") SR<-matrix(NA,nr=4,nc=4)#根据经济理论设定短期和长期结构系数矩阵 SR[4,2]<-0 LR<-matrix(NA,nr=4,nc=4) LR[1,2:4]<-0 LR[2:4,4]<-0 svec<-SVEC(vecm,LR=LR,SR=SR,r=1,lrtest=F,boot=T)#构建SVEC模型 summary(svec) SVEC EstimationResults: ========================  Call: SVEC(x = vecm, LR = LR, SR =SR, r = 1, lrtest = F, boot = T) Type:B-model  Sample size:81  Log Likelihood:-161.838  Number of iterations:12  Estimated contemporaneousimpact matrix:        prod       e        U    rw prod 0.58402  0.07434 -0.1525780.06900 e   -0.12029  0.26144 -0.1550960.08978 U    0.02526 -0.26720  0.0054880.04982 rw   0.11170  0.00000 0.483771 0.48791 Estimated standard errorsfor impact matrix:        prod      e      U    rw prod 0.09387 0.12013 0.216520.09304 e   0.06461 0.05874 0.16346 0.04542 U   0.05479 0.05526 0.06906 0.03374 rw  0.14126 0.00000 0.65297 0.09234 Estimated long run impactmatrix:        prod      e      Urw prod 0.7910  0.0000 0.0000  0 e    0.2024  0.5769 -0.4923 0 U   -0.1592 -0.3409  0.1408 0 rw  -0.1535  0.5961 -0.2495 0 Estimated standard errorsfor long-run matrix:       prod     e     Urw prod 0.1435 0.0000 0.0000 0 e   0.2220 0.1976 0.5540 0 U   0.1096 0.1032 0.1435 0 rw   0.17580.1806 0.2544  0 Covariance matrix of reducedform residuals (*100):        prod     e      U    rw prod 37.4642 -2.096 -0.2512 2.509 e   -2.0960 11.494 -6.9273 -4.467 U   -0.2512 -6.927  7.4544 2.978 rw   2.5087 -4.467  2.978348.457 LR[3,3]<-0 svec.oi<-update(svec,LR=LR,SR=SR,lrtest=T,boot=F)#模型的优化 svec.oi SVEC EstimationResults: ========================  Estimated contemporaneousimpact matrix:        prod      e       U     rw prod 0.58402  0.1199 0.11200 0.07039 e   -0.12029  0.3111 0.04089 0.09159 U    0.02526 -0.2571  0.096740.05083 rw   0.11170  0.0000 -0.473610.49778 Estimated long run impactmatrix:        prod      e      U rw prod 0.7910  0.0000 0.000000 0 e    0.2024  0.6768 0.238920 0 U   -0.1592 -0.3688 0.000000 0 rw  -0.1535  0.6456 0.003105 0 summary(svec.oi) SVEC EstimationResults: ========================  Call: SVEC(x = vecm, LR = LR, SR =SR, r = 1, lrtest = T, boot = F) Type:B-model  Sample size:81  Log Likelihood:-164.876  Number of iterations:17  LR overidentificationtest: LRoveridentification data: vecm Chi^2 = 6.1, df = 1, p-value= 0.01 Estimated contemporaneousimpact matrix:        prod      e       U     rw prod 0.58402  0.1199 0.11200 0.07039 e   -0.12029  0.3111 0.04089 0.09159 U    0.02526 -0.2571  0.096740.05083 rw   0.11170  0.0000 -0.473610.49778 Estimated long run impactmatrix:        prod      e      U rw prod 0.7910  0.0000 0.000000 0 e    0.2024  0.6768 0.238920 0 U   -0.1592 -0.3688 0.000000 0 rw  -0.1535  0.6456 0.003105 0 Covariance matrix of reducedform residuals (*100):        prod     e      U    rw prod 37.2953 -2.192 -0.1664 4.723 e   -2.1921 12.132 -7.4408 1.279 U   -0.1664 -7.441  7.8671-1.770 rw   4.7233  1.279 -1.769748.457 svec.irf<-irf(svec,response="U",n.ahead=48,boot=T) plot(svec.irf) fevd.U<-fevd(svec,n.ahead=48)$U matplot(fevd.U,type="l",lty=1:4)#失业率U的方差分解