/* Data 'a' Generated Process */
   proc iml;
      sig = {1.0  0.5, 0.5  1.25};
      phi = {1.2 -0.5, 0.6  0.3};
      call varmasim(y,phi) sigma = sig n = 100 seed = 46859;
      cn = {'y1' 'y2'};
      create a from y[colname=cn];
      append from y;
   run;;

   /* when the series has a linear trend */
   proc varmax data=a;
      model y1 y2 / p=1 trend=linear;
   run;
   
   /* Fit subset of AR order 1 and 3 */
   proc varmax data=a;
      model y1 y2 / p=(1,3);
   run;
   
   /* Check if the series is nonstationary */
   proc varmax data=a;
      model y1 y2 / p=1 dftest print=(roots);
   run;
   
   /* Fit VAR(1) in differencing */
   proc varmax data=a;
      model y1 y2 / p=1 print=(roots) dify=(1);
   run;
   
   /* Fit VAR(1) in seasonal differencing */
   proc varmax data=a;
      model y1 y2 / p=1 dify=(4) lagmax=5;
   run;
   
   /* Fit VAR(1) in both regular and seasonal differencing */
   proc varmax data=a;
      model y1 y2 / p=1 dify=(1,4) lagmax=5;
   run;
   
   /* Fit VAR(1) in different differencing */
   proc varmax data=a;
      model y1 y2 / p=1 dif=(y1(1,4) y2(1)) lagmax=5;
   run;
   
   /* Options related to prediction */
   proc varmax data=a;
      model y1 y2 / p=1 lagmax=3
                    print=(impulse covpe(5) decompose(5));
   run;
   
   /* Options related to tentative order selection */
   proc varmax data=a;
      model y1 y2 / p=1 lagmax=5 minic
                    print=(parcoef pcancorr pcorr);
   run;
   
   /* Automatic selection of the AR order */
   proc varmax data=a;
      model y1 y2 / minic=(type=aic p=5);
   run;
   
   /* Compare results of LS and Yule-Walker Estimators */
   proc varmax data=a;
      model y1 y2 / p=1 print=(yw);
   run;
   
   /* BVAR(1) of the nonstationary series y1 and y2 */
   proc varmax data=a;
      model y1 y2 / p=1
         prior=(lambda=1 theta=0.2 ivar);
   run;
   
   /* BVAR(1) of the nonstationary series y1 */
   proc varmax data=a;
      model y1 y2 / p=1
         prior=(lambda=0.1 theta=0.15 ivar=(y1));
   run;

   /* Data 'b' Generated Process */
   proc iml;
      sig = { 0.5  0.14 -0.08 -0.03,  0.14 0.71 0.16 0.1,
             -0.08 0.16  0.65  0.23, -0.03 0.1  0.23 0.16};
      sig = sig * 0.0001;
      phi = {1.2 -0.5 0.  0.1,  0.6 0.3 -0.2  0.5,
             0.4  0. -0.2 0.1, -1.0 0.2  0.7 -0.2};
      call varmasim(y,phi) sigma = sig n = 100 seed = 32567;
      cn = {'y1' 'y2' 'y3' 'y4'};
      create b from y[colname=cn];
      append from y;
   quit;

   /* Cointegration Rank Test using Trace statistics */
   proc varmax data=b;
      model y1-y4 / p=2 lagmax=4 cointtest;
   run;
   
   /* Cointegration Rank Test using Max statistics */
   proc varmax data=b;
      model y1-y4 / p=2 lagmax=4 cointtest=(johansen=(type=max));
   run;
   
   /* Common Trends Test using Filter(Differencing) statistics */
   proc varmax data=b;
      model y1-y4 / p=2 lagmax=4 cointtest=(sw);
   run;
   
   /* Common Trends Test using Filter(Residual) statistics */
   proc varmax data=b;
      model y1-y4 / p=2 lagmax=4 cointtest=(sw=(type=filtres lag=1));
   run;
   
   /* Common Trends Test using Kernel statistics */
   proc varmax data=b;
      model y1-y4 / p=2 lagmax=4 cointtest=(sw=(type=kernel lag=1));
   run;
   
   /* Cointegration Rank Test for I(2) */
   proc varmax data=b;
      model y1-y4 / p=2 lagmax=4 cointtest=(johansen=(iorder=2));
   run;
   
   /* Fit VECM(2) with rank=3 */
   proc varmax data=b;
      model y1-y4 / p=2 lagmax=4 print=(roots iarr)
                    ecm=(rank=3 normalize=y1);
   run;

   /* Weak Exogenous Testing for each variable */
   proc varmax data=b outstat=bbb;
      model y1-y4 / p=2 lagmax=4
                    ecm=(rank=3 normalize=y1);
      cointeg rank=3 exogeneity;
   run;
   
   /* Hypotheses Testing for long-run and adjustment parameter */
   proc varmax data=b outstat=bbb;
      model y1-y4 / p=2 lagmax=4
                    ecm=(rank=3 normalize=y1);
      cointeg rank=3 normalize=y1
         h=(1 0 0, 0 1 0, -1 0 0, 0 0 1)
         j=(1 0 0, 0 1 0, 0 0 1, 0 0 0);
   run;
   
   /* ordinary regression model */
   proc varmax data=grunfeld;
      model y1 y2 = x1-x3;
   run;
   
   /* Ordinary regression model with subset lagged terms */
   proc varmax data=grunfeld;
      model y1 y2 = x1 / xlag=(1,3);
   run;
   
   /* VARX(1,1) with no current time Exogenous Variables */
   proc varmax data=grunfeld;
      model y1 y2 = x1 / p=1 xlag=1 nocurrentx;
   run;
   
   /* VARX(1,1) with different Exogenous Variables */
   proc varmax data=grunfeld;
      model y1 = x3, y2 = x1 x2 / p=1 xlag=1;
   run;
   
   /* VARX(1,2) in difference with current Exogenous Variables */
   proc varmax data=grunfeld;
      model y1 y2 = x1 / p=1 xlag=2 difx=(1) dify=(1);
   run;