/*********     VAR Routines (Rapach&Wharff,June96)      ************************
*******************************************************************************
Author: David Rapach
==> For public, non-commercial use only
==> Please include proper attribution should this code be used in a
    research project or in other code
==> Author makes no performance guarantees
==> Program written for GAUSS Version 3.1.13
    (Vendor: Aptech Systems, Maple Valley, WA)
*******************************************************************************
******************************************************************************/

/******************************************************************************
Proc            : LR_LAG
Author          : David Rapach
Last revised    : 27 May 1996
Input           : z -- data matrix (z_1~z_2~etc, where z_i is column vector
                       of obs on ith endogenous variable)
                  p -- max lag (integer ge 2)
                  q -- significance level
Output          : Suggests VAR lag based on sequential testing
References      : Helmut Lutkepohl, _Introduction to Multiple Time Series
                      Analysis_ (New York: Springer-Verlag, 1993)
                  Christopher Sims, `Macroeconomics and reality,'
                      _Econometrica_ (Jan. 1980): 1-49
Notes           : -- Requires proc `varlags()' (provided)
                  -- Significance of each lag is tested (based on max lag
                     and significance level specified ) using the Sims
                     (1980) modified likelihood ratio statistic
                  -- If no lag ge 2 is significant, VAR(1) suggested
******************************************************************************/
proc(0)=lr_lag(z,p,q);
    local M,y,ylags,T,X,i,j,XL,bL,UL,sigL,XS,bS,US,sigS,c,LR,SLR,df,
          oldnv,pval;
    M=cols(z);                              @ # of endog vars @
    {y,ylags}=varlags(z,p);                 @ generating lags @
    T=rows(y);                              @ # of usable obs @
    i=0;
    j=0;
    oldnv=__fmtnv;
    __fmtnv="*.*lg"~2~2;
    "Sequential testing for VAR lag";
    do until ((i gt p-2) or (j ge 1));
        XL=ones(T,1)~ylags[.,1:p*M-i*M];    @ unrestricted data matrix @
        bL=y/XL;                            @ OLS estimates @
        UL=y-XL*bL;                         @ OLS resids @
        sigL=UL'UL/T;                       @ resids var-cov matrix @
        XS=ones(T,1)~ylags[.,1:p*M-((i+1)*M)];  @ restricted data matrix @
        bS=y/XS;                            @ OLS estimates @
        US=y-XS*bS;                         @ resids @
        sigS=US'US/T;                       @ resids var-cov matrix @
        c=cols(XL);                         @ Sims correction @
        LR=T*(ln(det(sigS))-ln(det(sigL))); @ LR stat @
        SLR=(T-c)*(ln(det(sigS))-ln(det(sigL)));@ modified LR stat @
        df=M*(cols(XL)-cols(XS));           @ # of restrictions @
        format 6,3;
        "--------------------------------------------------";
        "H0: VAR("$+ftocv(p-(i+1),1,0)$+")";
        "H1: VAR("$+ftocv(p-i,1,0)$+")";
        "------------";
        format 8,4;
        "LR stat         " LR  "P-value  " cdfchic(LR,df);
        pval=cdfchic(SLR,df);
        "Modified LR stat" SLR "P-value  " pval;
        if pval le q;
            j=1;
        else;
            i=i+1;
        endif;
    endo;
    "--------------------------------------------------";
    "Max lag            ";;call printfmt(p,1);?;
    "Significance level ";;call printfmt(q,1);?;
    "Sequential testing using the Sims (1980) modified";
    "likelihood-ratio statistic suggests a VAR("$+ftocv(p-i,1,0)$+")";
    __fmtnv=oldnv;
endp;

/****************************************************************************
Proc        : BOOTSTRP
Authors     : David Rapach and Jeffrey Wharff
Last revised: 6 January 1997           Previous Archived Version: June 1996
Input       : d     -- VAR data matrix (d_1~d_2~etc, where d_i is
                       column vector of observations on ith endogenous
                       variable)
              p     -- VAR order
              k     -- # of start-up transient observations
Output      : ystar -- VAR pseudo-data matrix (d_1star~d_2star~etc)
References  : * D.E. Runkle, "Vector autoregressions and reality,"
                    _Journal of Business and Economic Statistics_
                    (Oct. 1987): 437-42
              * J. Berkowitz and L. Kilian, "Recent developments
                in bootstrapping time series," Federal Reserve
                Board Finance and Economics Discussion Series
                Paper 1996-45 (November 1996)
Notes       : * Requires proc `varlags()' which is included
              * VAR pseudo-disturbances are generated by randomly
                    drawing from original disturbances
              * Intercepts included in VAR
              * Start-up transient observations are "extra" observations
                which allow initial observations to be random
                (see Berkowitz and Kilian 1996, p. 5)
****************************************************************************/
proc(1)=bootstrp(d,p,k);
    local M,y,X,T,bvar,U,segment,iter,z,q,e,N,Ustar,c,w,ystar,r,i,j;

    /* OLS VAR estimation */

    M=cols(d);                      @ # of endogenous variables @
    {y,X}=varlags(d,p);             @ generating lags @
    T=rows(y);                      @ # of usable observations @
    XX=ones(T,1)~X;                 @ adding intercepts @
    bvar=y/XX;                      @ OLS VAR parameters @
    U=y-XX*bvar;                    @ OLS residuals @

    /* Generating pseudo-disturbances */

    segment=zeros(T,1);             @ interval vector @
    iter=1;                         @ initializing @
    do until iter>T;                @ begin do loop @
        segment[iter,1]=iter*(1/T); @ creating intervals along (0,1] @
        iter=iter+1;                @ new iter @
    endo;                           @ end do loop @
    z=zeros(T+k,1);                 @ integer draw vector @
    q=1;                            @ initializing @
    do until q>T+k;                 @ begin do loop @
        e=rndu(1,1);                @ draw uniform random # @
        N=1;                        @ initializing @
        do until z[q,1]>0;          @ begin do loop @
            if e<=segment[N,1];     @ assigning integer to uniform draw @
            z[q,1]=N;               @ if draw <= Nth segment value, select @
            else;                   @ otherwise, go to next segment value @
            N=N+1;                  @ new N @
            endif;                  @ end if @
        endo;                       @ end do loop @
        q=q+1;                      @ new q @
    endo;                           @ end do loop @
    Ustar=zeros(T+k,M);             @ psuedo-disturbances @
    c=1;                            @ initializing @
    do until c>T+k;                 @ begin do loop @
        w=z[c,1];                   @ integer for cth obs @
        Ustar[c,.]=U[w,.];          @ pseudo-disturbance for cth obs @
        c=c+1;                      @ new c @
    endo;                           @ end do loop @

    /* Building pseudo-sample */

    ystar=zeros(T+k,M);                     @ pseudo realizations @
    r=X[1,.];                               @ lags for first observations @
    i=1;                                    @ initializing @
    do until i>T+k;                         @ begin do loop @
        ystar[i,.]=(1~r)*bvar+Ustar[i,.];   @ generating obs @
        r=shiftr(r,M,0);                    @ adjusting lags @
        j=1;                                @ initializing @
        do until j>M;                       @ begin do loop @
            r[.,j]=ystar[i,j];              @ filling in first lag @
            j=j+1;                          @ new j @
        endo;                               @ end do loop @
        i=i+1;                              @ new i @
    endo;                                   @ end do loop @
    ystar=d[1:p,.]|ystar;                   @ pseudo-sample @
    ystar=ystar[k+1:T+p+k,.];               @ adjusting sample size @
    retp(ystar);
endp;

/**********************  PROC VARLAGS  *****************************
**   Author: Alan G. Isaac
**   last update: 5 Dec 95      previous: 15 June 94
**   FORMAT
**          { x,xlags } = varlags(var,lags)
**   INPUT
**        var  - T x K matrix
**        lags - scalar, number of lags of var (a positive integer)
**   OUTPUT
**          x -     (T - lags) x K matrix, the last T-lags rows of var
**          xlags - (T - lags) x lags*cols(var) matrix,
**                  being the 1st through lags-th
**                  values of var corresponding to the values in x
**                  i.e, the appropriate rows of x(-1)~x(-2)~etc.
**   GLOBAL VARIABLES: none
**********************************************************************/
proc(2)=varlags(var,lags);
    local xlags;
    xlags = shiftr((ones(1,lags) .*. var)',seqa(1-lags,1,lags)
                                            .*. ones(cols(var),1),miss(0,0))';
    retp(trimr(var,lags,0),trimr(xlags,0,lags));
endp;
