THE FOLLOWING IS A QUICK GUIDE TO THE RASCHE/ROTEMBURG PROGRAMS AND
DATA FILES.  FOR A MORE DETAILED DESCRIPTION, PLEASE CHECK THE
INDIVIDUAL FILES.

FOR ROTEMBURG:

    JULIO.PGM:      HAS SAS PROGRAMS WHICH DO THE 4 REGRESSIONS.

    DATA.DAT:       CONTAINS DATA USED IN REGRESSIONS 1 AND 2.
                    THIS FILE IS ALSO AN ASCII TRANSFORMATION OF THE
                    DATA SET USED IN THE RASCHE PROGRAMS.

    JULIOQ.DAT:     CONTAINS DATA USED IN REGRESSIONS 3 AND 4.

    VELDATA.DAT:    CONTAINS DATA FOR FIGURE 2.  YOU CAN GET FIGURE 1
                    DATA FROM THE DATA.DAT FILE.

FOR RASCHE:

    M1VAR.PGM:      PROGRAM WHICH CALCULATES TABLES 1 AND 2 AS WELL AS
                    FIGURES 5-8.

    JOHANUM1.PGM:   PROGRAM WHICH CALCULATES TABLES 3 AND 4.

    DATA.DAT:       CONTAINS DATA USED IN THE RASCHE PAPER.  SINCE
                    THIS IS IN ASCII FORMAT, YOU WILL HAVE TO CHANGE
                    THE DATA STEPS IN THE RASCHE PROGRAMS TO READ IN
                    ASCII FORMATTED DATA, NOT RATS FORMATTED DATA.

    STLOUIS.APP:    PROGRAM WHICH CREATES THE CHARTS WHICH APPEAR IN THE
                    APPENDIX.

    NEWMON.RAT:     DATA FILE USED FOR THE PAPER.  THIS IS IN RATS
                    FORMAT.  THE DATA.DAT FILE IS AN ASCII TRANSFORM-
                    ATION OF THIS FILE.

    IMPULSE.UNR:    THE DATA WHICH WAS SAVED FROM A MONTE CARLO
                    SIMULATION AND PLOTTED IN THE APPENDIX.

************************************************************************

bma global 300 compile 3000
cal 55 1 4
all 0 91,4
ieval ibeg = cal(55,2)
ieval iend = (90:3)
ieval k = 4
ieval ibeg1 = ibeg+1
ieval ibegk = ibeg + k
*
*   THIS PROGRAM USES RATS VERSION 3.10 .
*
*   THIS PROGRAM IS USED TO PRODUCE TABLES 1 AND 2.  AS THIS PROGRAM
*   IS NOW WRITTEN, IT WILL PERFORM THE REGRESSION ANALYSIS OVER THE
*   1955.2 TO 1990.3 SAMPLE PERIOD.  TO RUN THE ANALYSIS OVER THE
*   1955.2 TO 1981.4 SAMPLE PERIOD, CHANGE THE IEND VARIABLE TO 81:4
*   AND ELIMINATE THE 1982 DUMMY VARIABLE, D82, FROM THE REGRESSION.
*
*   FIGURES 5-8 CAN ALSO BE CREATED IN THIS PROGRAM FROM THE IMPULSE
*   RESPONSE FUNCTION AT THE END OF THE PROGRAM.  THE PERTINENT FIGURES
*   ARE WHEN II = 3.
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί                 Data input and Transformation              Ί
*Ί           This section of the Program is not General       Ί
*Ί          Data should be read and transformed for your      Ί
*Ί                        Specific Problem                    Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
smpl ibeg iend
open data a:newmon.rat
data(format=rats) / m1q rtbq gnp gnp82
set d67 / = 0.0
set d67 67:4 91:4 = 1.0
set d79 / = 0.0
set d79 79:4 91:4 = 1.0
set d82  / = 0.0
set d82 82,1 91:4 = 1.0
set skip / = 1.0
smpl ibeg iend
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί     Construct levels of Variables here if Transformations  Ί
*Ί                   are Required (p variables)               Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
set qp / =  log(gnp(t))- log(gnp82(t))
set x3 / = log(gnp82(t))
set x2 / =  400*(qp(t) -qp(t-1))
set x1 / =  log(m1q(t))
set x4 / = rtbq(t)
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί          End of Data Input and Transformation Section       Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
smpl ibegk iend
system 1 to 4
variables x3 x2 x1 x4
lags 1 to k
det constant d67 d79 d82
end(system)
estimate(outsigma=omega)
declare vect std(4)
do ii = 1,4
   eval std(ii) = 1.0/sqrt(omega(ii,ii))
   end do ii
matrix omega = diag(std)*omega*diag(std)
write omega
declare rect remp(4,4)
matrix remp = decomp(omega)
set imp1 1 41 = 0.0
set imp2 1 41 = 0.0
set imp3 1 41 = 0.0
set imp4 1 41 = 0.0
set imp5 1 41 = 0.0
set imp6 1 41 = 0.0
set imp7 1 41 = 0.0
set imp8 1 41 = 0.0
do ii = 1,4
   IMPULSE(DECOMP=REMP,NOPRINT) 4 40 II
   #1 imp1 2
   #2 imp2 2
   #3 imp3 2
   #4 imp4 2
   * Price level
   set imp8 2 41 = .0025*imp2(t) + imp8(t-1)
   * nominal balances
   set imp5 2 41 = imp3(t) + imp8(t)
   * velocity growth
   set imp6 2 41 =  400*(imp1(t) - imp1(t-1)) + imp2(t) - $
   400*(imp5(t) - imp5(t-1))
   * nominal income growth
   set imp7 2 41 = 400*(imp1(t) - imp1(t-1)) + imp2(t)
   print(nodates) 2 41 imp1 imp2 imp3 imp4 imp5 imp6 imp7 imp8
   end do ii
end

************************************************************************

bma global 300 compile 3000
output noregress noqstat
*
*   THIS PROGRAM USES RATS VERSION 3.10 .
*
*   THIS PROGRAM CAN BE USED TO CREATE TABLES 3 AND 4.  EVERYTHING
*   THAT ONE NEEDS IS ALREADY LABELLED THROUGHOUT THE PROGRAM.  NONE
*   OF THE FIGURES WAS CREATED WITH THIS PROGRAM.  ALL THAT NEEDS TO
*   BE DONE IS TO CHANGE THE SAMPLE PERIODS, AS SPECIFIED BY THE IBEG
*   AND IEND VARIABLES, TO LOOK AT A SPECIFIC PERIOD.  ALSO, TWO X2
*   VARIABLES ARE CREATED, SO REMEMBER TO COMMENT OUT ONE OF THEM
*   DEPENDING ON WHICH TABLE YOU ARE RE-CREATING.  ALSO, WHEN TESTING
*   EARLIER SAMPLE PERIODS, YOU HAVE TO ELIMINATE SOME OF THE SAMPLE
*   DEPENDENT DUMMY VARIABLES, IE D82 AND D79.
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί    Maximum Likelihood Estimation of Vector Error Correction Ί
*Ί                            Models                           Ί
*Ί                                                             Ί
*Ί              Estimation Procedure of S. Johansen            Ί
*Ί         Journal of Economic Dynamics and Control, 1988      Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
cal 55 1 4
all 0 91,4
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί                Set Parameters of Johansen Model             Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
ieval p = 4
ieval s = 3
ieval k = 4
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί      set number of cointegration vectors to compute alpha   Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
ieval r = 2
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί    Set number of columns for Constraint Matrix (A) on alpha Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
ieval mr = p
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί                    Set p1 = p+1 to test H*2                 Ί
*Ί            Johansen and Juselius (1989) Section II          Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
ieval p1 = p
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί               Set Sample Period for Regressions             Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
ieval ibeg = cal(55,2)
ieval iend = (90:3)
*
if mr < p
    ieval mr1 = p-mr
  else
    ieval mr1 = p
  end if
declare vect test(s)
declare symm s00(p,p) skk(p1,p1) m(s,s) mm(s,s) sd1(p-r,p-r)
declare rect s0k(p,p1) c(s,s) seivec(s,s) ci(s,s) h(p1,s) tri(s,s)
declare symm skkb(p1,p1) saab(mr,mr) b00b(mr1,mr1)
declare rect aj(p,mr) bj(p,mr1) sakb(mr,p1) g0 b0(r,1) mu0(p,4)
declare rect bet gam1 mu aprp(p,p-r) one(s,1) pi11(p,p) b0pr(p,4)
declare rect psi(p,4) betaprp a1 d1 v02 vk2 mubar(p-r,4)
ewise tri(i,j) = (i<=j)*(-1.0)
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί     Set H matrix for restrictions on Cointegration Vector   Ί
*Ί            Johansen and Juselius (1989) Section V           Ί
*Ί                   Default = Identity Matrix                 Ί
*Ί                                                             Ί
*Ί                 Example of Restriction Matrix               Ί
*Ί                           (p=4,s=2)                         Ί
*Ίinput h                                                      Ί
*Ί 1 0                                                         Ί
*Ί 0 1                                                         Ί
*Ί-1 0                                                         Ί
*Ί 0 0                                                         Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
*matrix h = iden(1.0)
input h
 1 0 0
 0 1 0
-1 0 0
 0 0 1
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί     Set A and B matricies for restrictions on alpha matrix  Ί
*Ί            Johansen and Juselius (1989) Section VI          Ί
*Ί                                                             Ί
*Ί                Example of Restriction Matricies             Ί
*Ί                           (p=4,mr=3)                        Ί
*Ί input aj                                                    Ί
*Ί 1 0 0                                                       Ί
*Ί 0 1 0                                                       Ί
*Ί 0 0 1                                                       Ί
*Ί 0 0 0                                                       Ί
*Ί input bj                                                    Ί
*Ί 0                                                           Ί
*Ί 0                                                           Ί
*Ί 0                                                           Ί
*Ί 1                                                           Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
matrix aj = iden(1.0)
matrix bj = iden(0.0)
ieval k1 = k-1
ieval ibeg1 = ibeg+1
ieval ibegk = ibeg+k
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί                 Data input and Transformation              Ί
*Ί           This section of the Program is not General       Ί
*Ί          Data should be read and transformed for your      Ί
*Ί                        Specific Problem                    Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
smpl ibeg iend
open data a:newmon.rat
data(format=rats) / m1q rtbq gnp gnp82
set d67 / = 0.0
set d67 67:4 91:4 = 1.0
set d79 / = 0.0
set d79 79:4 91:4 = 1.0
set d82 / = 0.0
set d82 82,1 91:4 = 1.0
set skip / = 1.0
smpl ibeg iend
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί     Construct levels of Variables here if Transformations  Ί
*Ί                   are Required (p variables)               Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
set qp / =  log(gnp(t))- log(gnp82(t))
set x1 / = log(m1q(t)) -qp(t)
*set x2 / =  400*(qp(t) -qp(t-1))
set x2 / = 400.*(log(m1q(t))-log(m1q(t-1)))
set x3 / = log(gnp82(t))
*set x4 / = log(rtbq(t))
set x4 / = rtbq(t)
smpl ibeg1  iend
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί         Construct First Differences of Variables Here      Ί
*Ί                         (p variables)                      Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
set dx1 / = x1(t) - x1(t-1)
set dx2 / = x2(t) - x2(t-1)
set dx3 / = x3(t) - x3(t-1)
set dx4 / = x4(t) - x4(t-1)
smpl ibegk iend
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί        Construct Lagged (t-k) Levels of Variables Here     Ί
*Ί                         (p variables)                      Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
set x1k / = x1(t-k)
set x2k / = x2(t-k)
set x3k / = x3(t-k)
set x4k / = x4(t-k)
ieval ibegk = ibegk+1
smpl ibegk iend
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί          End of Data Input and Transformation Section       Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
*
*
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί                   VAR Regression Section II                 Ί
*Ί            Regressions are specific to your problem         Ί
*Ί  All Regression Residuals and Coefficients must be RetrievedΊ
*Ί                        (p1 Regressions)                     Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
linreg(smpl=skip)  x1k / res1 coef1
#constant d67 d79 d82 -dx1 1 k1 -dx2 1 k1 -dx3 1 k1 -dx4 1 k1
eval sqt = sqrt(float(nobs))
eval tt = float(nobs)
eval ttinv = 1.0/tt
ieval pk1 = nreg
linreg(smpl=skip)  x2k / res2 coef2
#constant d67 d79 d82 -dx1 1 k1 -dx2 1 k1 -dx3 1 k1 -dx4 1 k1
linreg(smpl=skip)  x3k / res3 coef3
#constant d67 d79 d82 -dx1 1 k1 -dx2 1 k1 -dx3 1 k1 -dx4 1 k1
linreg(smpl=skip)  x4k / res4 coef4
#constant d67 d79 d82 -dx1 1 k1 -dx2 1 k1 -dx3 1 k1 -dx4 1 k1
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί   Make Matricies from Regression Residuals and Coefficients Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
make resk
#res1 res2 res3 res4
make(transpose) vk 1 pk1
#coef1 coef2 coef3 coef4
overlay vk(1,2) with vk2(p,1)
clear res1 resq2 res3 coef1 coef2 coef3
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί                   VAR Regression Section I                  Ί
*Ί            Regressions are specific to your problem         Ί
*Ί  All Regression Residuals and Coefficients must be RetrievedΊ
*Ί                         (p Regressions)                     Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
linreg(smpl=skip)  dx1 / res1 coef1
#constant d67 d79 d82 -dx1 1 k1 -dx2 1 k1 -dx3 1 k1 -dx4 1 k1
linreg(smpl=skip) dx2 / res2  coef2
#constant d67 d79 d82 -dx1 1 k1 -dx2 1 k1 -dx3 1 k1 -dx4 1 k1
linreg(smpl=skip) dx3 / res3  coef3
#constant d67 d79 d82 -dx1 1 k1 -dx2 1 k1 -dx3 1 k1 -dx4 1 k1
linreg(smpl=skip) dx4 / res4  coef4
#constant d67 d79 d82 -dx1 1 k1 -dx2 1 k1 -dx3 1 k1 -dx4 1 k1
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί   Make Matricies from Regression Residuals and Coefficients Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
declare rect betpr(r,p) alfa(p,r) bigpi(p,p)
declare rect eye(s,s)  jj(p*p,p*(p+1)/2)
declare rect rh1(1,2) rh2(1,1)
declare rect gamma(p,pk1)
declare symn m1rr(p,p) sigma(p*(p+1)/2,p*(p+1)/2) omega(p,p) biglm(p,p)
make res0
#res1 res2 res3 res4
make(transpose) v0 1 pk1
#coef1 coef2 coef3 coef4
overlay v0(1,2) with v02(p,1)
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί       Estimation of Cointegration Vector and Hypothesis     Ί
*Ί                            Testing                          Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
output noecho
matrix s00 = tr(res0)*res0
if mr < p
    matrix b00b = inv(tr(bj)*s00*bj)
  else
    matrix b00b = iden(1.0)
  end if
matrix skk = tr(resk)*resk
matrix s0k = tr(res0)*resk
matrix skkb = skk - tr(s0k)*bj*b00b*tr(bj)*s0k
matrix saab = tr(aj)*s00*aj - tr(aj)*s00*bj*b00b*tr(bj)*s00*aj
matrix sakb = tr(aj)*s0k - tr(aj)*s00*bj*b00b*tr(bj)*s0k
matrix mm = tr(h)*skkb*h
matrix c = decomp(mm)
matrix m = inv(c)*tr(h)*tr(sakb)*inv(saab)*sakb*h*inv(tr(c))
display '                  Eigenvalues'
display ' '
eigen m eival eivec
write eival
matrix seivec = sqt*inv(tr(c))*eivec
display '             Eigenvectors (columnwise)'
display ' '
write seivec
display '           Number of Observations (T) = ' tt
display ' '
display '           Johansen Trace Test Statistics'
display 'Ho:   r=0          r<=1          r<=2         r<=3'
display ' '
dimension one(s,1)
ewise one(i,j) = 1.0
matrix test = tt*tri*log(one-eival)
write test
matrix eye = iden(1.0)
overlay seivec(1,1) with bet(s,r)
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί     extract beta' [rxs] matrix from matrix of eigenvectors  Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
matrix betpr = tr(h*bet)
display '                 beta prime Matrix'
display ' '
write betpr
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί                      compute alpha [sxr]                    Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
matrix alfa = scale(-ttinv)*aj*inv(tr(aj)*aj)*sakb*tr(betpr)
display '                   alpha Matrix'
display ' '
write alfa
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί                    compute big lambda [pxp]                 Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
matrix biglm =ttinv*s00-alfa*tr(alfa)
display '                 Lambda Matrix'
display ' '
write biglm
matrix lndet = log(det(biglm))
display '     -2*(Log of Maximized Likelihood) = det(Lambda)'
display ' '
write lndet
input jj
1 0 0 0 0 0 0 0 0 0
0 1 0 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0 0
0 0 0 1 0 0 0 0 0 0
0 1 0 0 0 0 0 0 0 0
0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0
0 0 0 0 0 0 1 0 0 0
0 0 1 0 0 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0
0 0 0 0 0 0 0 1 0 0
0 0 0 0 0 0 0 0 1 0
0 0 0 1 0 0 0 0 0 0
0 0 0 0 0 0 1 0 0 0
0 0 0 0 0 0 0 0 1 0
0 0 0 0 0 0 0 0 0 1
matrix omega = inv(biglm)
matrix sigma = inv(scale(tt/2.0)*tr(jj)*kroneker(omega,omega)*jj)
display '       Variance-Covariance Matrix of Covariances'
display ' '
write sigma
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί                  compute big pi matrix  [pxp]               Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
matrix bigpi = alfa*betpr
display '                   Pi Matrix'
display ' '
write bigpi
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί   compute matrix of VAR regression coefficients [pxp(k-1)+1]Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
matrix gamma = v0+bigpi*vk
display '                  Gamma Matrix'
display ' '
write gamma
ieval parms = r*s
declare rect d(r,r) betac(s,r) w1(r,r) w2a(s,s) betacvr(parms,parms)
declare rect vbetac pltb(r,r) betacp(r,s) pull restri(parms,2) wald(1,1)
declare rect vbet(parms,1) betapcp pulh(s,r)
dimension eye(s,s)
matrix eye = iden(1.0)
overlay eye(1,1) with pull(s,r)
matrix pulh = pull
matrix d = inv(diag(tr(pull)*eival))
matrix pltb = tr(pull)*bet
matrix betac = bet*inv(pltb)
matrix betacp = tr(betac)
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί                    Decompose Constant Term                  Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
overlay gamma(1,1) with mu(p,4)
overlay gamma(1,pk1+1-p*(k-1)) with gam1(p,p*(k-1))
dimension m(p,p) eye(p,p) one(k-1,1)
ewise one(i,j) = 1.0
matrix eye = iden(1.0)
matrix m = eye - tr(betpr)*inv(betpr*tr(betpr))*betpr
eigen m meival meivec
overlay meivec(1,1) with betaprp(p,p-r)
matrix m1rr = inv(biglm)* $
    (eye-alfa*inv(tr(alfa)*inv(biglm)*alfa)*tr(alfa)*inv(biglm))
eigen m1rr mmval mmvec
overlay mmvec(1,1) with a1(p,p-r)
overlay mmval(1) with d1(p-r,1)
matrix sd1 = diag(d1)
matrix aprp = a1*decomp(sd1)
matrix pi11 = eye -gam1*kroneker(eye,one) - scale(k-2)*bigpi
matrix c = betaprp*inv(tr(aprp)*pi11*betaprp)*tr(aprp)
display '             C(1) Matrix'
display ' '
write c
matrix psi = c*mu
display '     Trends in Nonstationary Variables'
display ' '
write psi
matrix b0pr = inv(tr(alfa)*alfa)*tr(alfa)*(mu-(eye-gam1*kroneker(eye,one))*psi)
display '        Means of Cointegration Vectors'
display ' '
write b0pr
display '        Beta Prp -- KPSW A(o)'
display ' '
matrix m = eye - h*tr(betacp)*inv(betacp*tr(betacp))*betacp*tr(h)
eigen m meival meivec
overlay meivec(1,1) with betapcp(p,p-r)
declare rect betaptr(p-r,p) az
matrix betaptr = tr(betapcp)
overlay betaptr(1,r+1) with az(p-r,p-r)
matrix betapcp = tr(inv(az)*betaptr)
write betapcp
overlay betacp(1,1) with vbetac(parms,1)
dimension eye(r,r)
matrix eye = iden(1.0)
matrix w1 = inv((tr(pulh)*bet)*inv(d-eye)*(tr(bet)*pulh))
dimension eye(s,s)
matrix  eye = iden(1.0)
matrix w2a = (eye-betac*tr(pulh))*inv(mm)*(eye-pulh*tr(betac))
matrix betacvr = kroneker(w2a,w1)
display '              Normalized Cointegration Vectors'
display ' '
write betac
display '          Covariance Matrix of Normalized Cointegration Vector'
display ' '
write betacvr
display '              Normalized Means of Cointegration Vector'
display ' '
matrix b0pr = inv(tr(pltb))*b0pr
write b0pr
display '              Normalized Alpha Matrix'
display ' '
matrix alfa = alfa*tr(pltb)
write alfa
ieval nrest=1
declare rect krest(nrest,parms) chi(1,1) vbetar(parms,1) b0prr(2,r)
declare symn betacvrr(parms,parms)
matrix b0prr = tr(b0pr)*inv(pltb)
write b0prr
input krest
0 0 0 1 0 1
matrix chi = tr(krest*vbetac)*inv(krest*betacvr*tr(krest))*krest*vbetac
cdf chisq  chi(1,1)  nrest
matrix vbetar = vbetac - $
 betacvr*tr(krest)*inv(krest*betacvr*tr(krest))*krest*vbetac
display ' Constrained Cointegrating Vectors'
display ' '
write vbetar
matrix betacvrr = betacvr - $
 betacvr*tr(krest)*inv(krest*betacvr*tr(krest))*krest*betacvr
display ' Covariance Matrix of Constrained Cointegrating Vectors'
display ' '
write betacvrr
*
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί           Compute WALD tests of Linear Restrictions                     Ί
*Ί              On Cointegration Vectors if Desired                        Ί
*Ί                                                                         Ί
*Ί       Example  (nr = the Number of Linear Restrictions)                 Ί
*Ί                (kr = Matrix of Restriction Coefficients)                Ί
*Ίieval nr = 1                                                             Ί
*Ίdeclare rect kr(parms,nr)                                                Ί
*Ίinput kr                                                                 Ί
*Ί0                                                                        Ί
*Ί1                                                                        Ί
*Ί0                                                                        Ί
*Ί0                                                                        Ί
*Ίmatrix wald = tr(vbetac)*kr*inv(tr(kr)*betacvr*kr)*tr(kr)*vbetac         Ί
*Ίdisplay '     Wald Test for Linear Restrictions on Cointegration Vectors'Ί
*Ίcdf=chisqr wald nr                                                       Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
*
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί         compute matrix of regression residuals [ttxp]       Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
declare integer iresid
scratch(matrix=resid) p ibegk iend iresid
matrix resid = const(0.0)
matrix resid = res0 + resk*tr(bigpi)
display '            Matrix of Regression Residuals'
display ' '
write resid
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί     move residuals to series and compute autocorrelations   Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
do ii = 1,p
   statistics(higher) iresid+ii
   correlate iresid+ii / 12 corrs
   eval qstat = 0.0
   do ij = 1,12
      eval qstat=qstat+tt*(tt+2)*corrs(ij+1)**2/(tt-ij)
      end do ij
   cdf chisqr qstat 12
   end do ii
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί                     Supplementary Sections                              Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί    Compute Coefficients of a VAR in Levels for Forecasting  Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
 declare rect t2(p,pk1+p) t3(p,pk1+p) tl2 tl3 t1a
 declare rect t1p(pk1+p,pk1) t1p1 t11(p*(k-1),pk1+p)
 declare rect t1(p*(k-1)+4,pk1+p) matcoef(p,pk1+p) mcoefpr(pk1+p,p)
 declare vector over1 over2 over3 over4
 dimension eye(p,p)
 matrix eye = iden(1.0)
 ieval ps = pk1-p*k1+1
 overlay t2(1,ps) with tl2(p,p*k)
 overlay t3(1,ps) with tl3(p,p*k)
 overlay t11(1,ps) with t1a(p*(k-1),p*k)
 overlay t1p(1,ps) with t1p1(p*k+ps-1,p*(k-1))
 ewise t2(i,j) = 0.0
 ewise t3(i,j) = 0.0
 declare rect row(1,k)
 ewise row(i,j) = 0
 eval row(1,k) = 1.0
 matrix tl2 = kroneker(eye,row)
 ewise row(i,j) = 0.0
 eval row(1,1) = 1.0
 matrix tl3 = kroneker(eye,row)
 dimension row(k-1,k)
 ewise t11(i,j) = 0.0
 ewise row(i,j) = 0.0
 do ik = 1,k1
    eval row(ik,ik) = 1.0
    eval row(ik,ik+1) = -1.0
    end do ik
 matrix t1a = kroneker(eye,row)
 ewise t1p(i,j) = 0.0
 do ik = 1,ps
     eval t1p(ik,ik) = 1.0
     end do ik
 matrix t1p1 = tr(t11)
 matrix t1 = tr(t1p)
matrix matcoef = gamma*t1-bigpi*t2+t3
matrix mcoefpr = tr(matcoef)
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί       Define equations for Input to Forecasting Routines    Ί
*Ί                         (p equations)                       Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
equation 1 x1
#constant d67 d79 d82 x1{1 to k} x2{1 to k} x3{1 to k} x4{1 to k}
equation 2 x2
#constant d67 d79 d82 x1{1 to k} x2{1 to k} x3{1 to k} x4{1 to k}
equation 3 x3
#constant d67 d79 d82 x1{1 to k} x2{1 to k} x3{1 to k} x4{1 to k}
equation 4 x4
#constant d67 d79 d82 x1{1 to k} x2{1 to k} x3{1 to k} x4{1 to k}
  overlay mcoefpr(1,1) with over1(pk1+p)
  overlay mcoefpr(1,2) with over2(pk1+p)
  overlay mcoefpr(1,3) with over3(pk1+p)
  overlay mcoefpr(1,4) with over4(pk1+p)
  associate(perm) 1 over1
  associate(perm) 2 over2
  associate(perm) 3 over3
  associate(perm) 4 over4
display        ' Coefficient Vectors of VAR in Levels'
display ' '
write mcoefpr
declare symn bigcorr(p,p)
declare vect std(p)
do ii = 1,p
    eval std(ii) = 1.0/sqrt(biglm(ii,ii))
    end do ii
matrix bigcorr = diag(std)*biglm*diag(std)
do ii = 1,p
    impulse p 40 ii bigcorr
    #3
    #2
    #1
    #4
    end do ii
* Common Trends Analysis and Decompositions
declare rect f(p-r,p) r3(p-r,p) r1(p,p) phi(p-r,p-r)
declare rect norm1(p-r,p-r) ahat(p,p-r) fstar(p-r,p) r0(p,p)
declare rect hs(p,p) hsinv(p,p) hsinvu(p,p)
declare symn sigeta(p-r,p-r)
dimension h(p,p)
* Matrix to change columns in A(0) if necessary
input phi
0 1
1 0
matrix betapcp = betapcp*phi
display '        KPSW A(0) Matrix after column exchange'
display ' '
write betapcp
matrix f = inv(tr(betapcp)*betapcp)*tr(betapcp)*c
write f
matrix sigeta = f*biglm*tr(f)
matrix phi = decomp(sigeta)
matrix norm1 = iden(1.0)
ieval pr = p-r
do ik = 1,pr
    eval norm1(ik,ik) = 1.0/phi(ik,ik)
    end do ik
matrix phi = phi*norm1
display '        KPSW Pi Matrix'
display ' '
write phi
cntrl(pause)
matrix ahat = betapcp*phi
display '        KPSW Ahat Matrix'
display ' '
write ahat
matrix fstar = inv(tr(ahat)*ahat)*tr(ahat)*c
matrix mm = tr(c)*c
eigen mm eival hinv
matrix h = inv(hinv)
matrix r0 = c*hinv
matrix r3 = fstar*hinv
matrix r1 = iden(1.0)
do ii = 1, pr
do ij = 1, pr
    eval r1(ii,ij) = r3(ii,ij)
    end do ij
    end do ii
matrix hs = r1*h
matrix hsinv = inv(hs)
dimension sigeta(p,p)
matrix sigeta = hs*biglm*tr(hs)
matrix r0 = decomp(sigeta)
write r0
matrix r1 = iden(1.0)
do ii = 1, p
    eval r1(ii,ii) = r0(ii,ii)
    end do ii
matrix r0 = inv(r0)
matrix hs = r1*r0*hs
matrix hsinv = inv(hs)
display '         KPSW Gamma(0) Matrix'
display ' '
write hsinv
matrix r0 = c*hsinv
do ii = 1, p
    impulse(decomp=hsinv) p 40 ii
    #1
    #2
    #3
    #4
    end do ii
matrix sigeta = hs*biglm*tr(hs)
display '   Covariance Matrix of Sturctural Disturbances'
display '          Check: Should be Diagonal'
display ' '
write sigeta
matrix r0 = decomp(sigeta)
matrix hsinvu = hsinv*r0
errors(decomp=hsinvu) p 40
#1
#2
#3
#4
end
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί         Construct Simulations or Forecasts if Desired       Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
smpl 90:1  91:3
steps(print) 3 7 90:1
# 1 x1e
# 2 x2e
# 3 x3e
set qvele / = x2e(t) - x1e(t)
print / qvele
end
graph 2
#x1
#x1e
graph 2
#x2
#x2e
graph 2
#x3
#x3e
eval biglm(2,1) = 0.0
eval biglm(3,1) = 0.0
write biglm
end
matrix rh1 = tr(alfa)*(eye-gam1*kroneker(eye,one))*psi
matrix rh2 = tr(alfa)*bigpi*vk2
eval rh3 = rh1(1,2) - rh2(1,1)
display rh3
eval rest = (rh3-alfa(1,1)*v02(1,1)-alfa(2,1)*v02(2,1))/alfa(3,1)
display v02(3,1)
display rest
linreg(smpl=skip) dx3 / res3  coef3
#constant  -dx1 1 k1 -dx2 1 k1 -dx3 1 k1
restrict(create) 1 res3 coef3
#2
#1 rest
end do ik

************************************************************************

bma global 300 compile 3000
cal 55 1 4
all 0 91,4
ieval ibeg = cal(55,2)
ieval iend = (90:3)
ieval p = 4
ieval k = 4
ieval r = 2
ieval icon = 4
ieval nstep = 40
ieval nvar = 12
ieval ndraws = 1
ieval ibeg1 = ibeg+1
ieval ibegk = ibeg + k
ieval k1 = k-1
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*  THIS PROGRAM USES RATS VERSION 3.10 .
*
*  THIS PROGRAM IS USED TO RE-CREATE TABLE 5, FIGURES 9-22, AND THE
*  MONTE CARLO SIMULATIONS WHICH APPEAR IN THE APPENDIX.  TABLE 5
*  IS ALREADY SET UP.  FOR THE 1981 DATA FOR THE FIGURES, YOU HAVE
*  TO CHANGE A FEW THINGS.  ONE, CHANGE THE END OF THE SAMPLE PERIOD
*  TO 81:4.  TWO, THE CIV1 VARIABLE HAS TWO VALUES, ONE FOR THE 81:4
*  ENDING PERIOD AND ANOTHER FOR THE 90:3 ENDING PERIOD.  CIV1 WILL
*  HAVE THE .1201614 CO-INTEGRATING VECTOR IF THE SAMPLE ENDING PERIOD
*  IS 90:3.  OTHERWISE, USE THE CIV1 VARIABLE WITH THE .1217731
*  CO-INTEGRATING VECTOR.  THERE ARE TWO OTHER PLACES IN THE PROGRAM
*  WHERE THESE TWO VALUES ARE INTER-CHANGED, SO MAKE SURE THAT YOU
*  ARE AWARE OF THE SAMPLE ENDING PERIOD.  ALSO, AS BEFORE, DON'T
*  FORGET TO ELIMINATE THE D82 VARIABLE.  FURTHERMORE, YOU HAVE TO
*  CHANGE THE VALUE OF ICON TO 3, NOT 4.
*
*  YOU CAN ALSO PLOT VARIOUS
*  SERIES FOR SOME OF THE FIGURES.  RELEVANT SERIES INCLUDE THE
*  FOLLOWING:  #35 #36 #37 #40 #41 #47 #48 #49 #52 #53.  IF YOU CHANGE
*  THE NOPRINT OPTION TO PRINT IN THE IMPULSE COMMAND, YOU WILL GET A
*  48 X 40 MATRIX WHICH SHOULD INCLUDE THIS DATA.  IN ADDITION, YOU WILL
*  ALSO GET THE DATA FOR THE TRANSITORY SHOCKS WHICH ARE PLOTTED IN
*  FIGURES 19 AND 20.  FIGURES 21 AND 22 ARE CALCULATED USING THE
*  TRANSITORY SHOCK DATA AND THE OLS COEFFICIENTS AS SHOWN IN THE PAPER.
*  IF YOU WANT TO RUN THE MONTE CARLO SIMULATION, CHANGE THE VALUE OF
*  THE VARIABLE NDRAWS FROM 1 TO 1000 AND RUN.
*
*Ί                 Data input and Transformation              Ί
*Ί           This section of the Program is not General       Ί
*Ί          Data should be read and transformed for your      Ί
*Ί                        Specific Problem                    Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
smpl ibeg iend
open data a:newmon.rat
data(format=rats) / m1q rtbq gnp gnp82
set d67 / = 0.0
set d67 67:4 91:4 = 1.0
set d79 / = 0.0
set d79 79:4 91:4 = 1.0
set d82  / = 0.0
set d82 82,1 91:4 = 1.0
set skip / = 1.0
smpl ibeg iend
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί     Construct levels of Variables here if Transformations  Ί
*Ί                   are Required (p variables)               Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
set qp / =  log(gnp(t))- log(gnp82(t))
set x1 / = log(m1q(t)) -qp(t)
set x2 / = 400.*(log(m1q(t)) - log(m1q(t-1)))
set x3 / = log(gnp82(t))
set x4 / = rtbq(t)
set x5 / =   400*(qp(t) -qp(t-1))
set x6 / = x3(t) - x1(t)
set civ1 / = X1(T)-X3(T)+.1201614*X4(T)
*set civ1 / = X1(T)-X3(T)+.1217731*X4(T)
set civ2 / = x2(t)-x4(t)
smpl ibeg1 iend
set dx1 / = x1(t) - x1(t-1)
set dx2 / = x2(t) - x2(t-1)
set dx3 / = x3(t) - x3(t-1)
set dx4 / = x4(t) - x4(t-1)
smpl ibegk iend
*ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
*Ί          End of Data Input and Transformation Section       Ί
*ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
system 1 to 4
variables dx1 dx2 dx3 dx4
lags 1 to k1
kfset xxx
det constant d67 d79 civ1{4} civ2{4}
end(system)
declare integer icoef iresid
ieval pk1 = p*k1
ieval ncoef = pk1+icon+2
scratch(matrix=coeff) p 1 ncoef icoef
scratch(matrix=resid) p ibegk iend iresid
eqv iresid+1 to iresid+4
resid1 resid2 resid3 resid4
estimate(outsigma=biglm) ibegk iend  iresid+1 icoef+1
display '     Covariance Matrix of VECM Errors'
display ' '
write biglm
associate(nocoeffs) 1 0  resid1
associate(nocoeffs) 2 0  resid2
associate(nocoeffs) 3 0  resid3
associate(nocoeffs) 4 0  resid4
equation(identity) 5 x1
# x1{1} dx1
associate(perm) 5
# 1.0 1.0
equation(identity) 6 x2
#x2{1} dx2
associate(perm) 6
# 1.0 1.0
equation(identity) 7 x3
#x3{1} dx3
associate(perm) 7
# 1.0 1.0
equation(identity) 8 x4
#x4{1} dx4
associate(perm) 8
# 1.0 1.0
equation(identity) 9 civ1
# x1 x3 x4
associate(perm) 9
# 1.0 -1.0 .1201614
equation(identity) 10 civ2
# x2 x4
associate(perm) 10
# 1.0 -1.0
equation(identity) 11 x5
# x2 x1{0 to 1}
associate(perm) 11
# 1.0 -400.0 400.0
equation(identity) 12 x6
# x3 x1
associate(perm) 12
# 1.0 -1.0
declare rect sxx(ncoef,ncoef) sinv(p,p) swish(p,p) ranw(p,nobs) $
   rantr(nobs,p) ranc(ncoef,p) sxxuover  betap(2,p) gamm0 gamma  $
   bigphi alfa bigpi(p,p) alpha vcoefp(p,pk1+icon+2) gam(p,pk1)
declare rect first(nstep,nvar*p) second(nstep,nvar*p)
declare vect betadraw betaols(ncoef*p) sxxu(ncoef*p)
declare rect remp(p,p)
declare vect std(p) mmval(p)
declare rect betac(p,r) pltb(r,r) betacp(r,p) pull betapcp pulh(p,r) h(p,p)
declare rect eye(p,p) mmvec(p,p) meivec(p,p)
declare symn m(p,p) sd1(p-r,p-r) m1rr(p,p) vmat(p,p)
declare rect betaprp a1 d1 aprp(p,p-r) pi11(p,p) c(p,p) one(k-1,1)
declare rect f(p-r,p) r3(p-r,p) r1(p,p) phi(p-r,p-r) hld(p-r,p-r)
declare rect norm1(p-r,p-r) ahat(p,p-r) fstar(p-r,p) r0(p,p)
declare rect hs(p,p) hsinv(p,p) hsinvu(p,p) kpswpim(p-r,p-r) kpswpiv(p-r,p-r)
declare symn sigeta(p-r,p-r) mm(p,p)
declare rect betaptr(p-r,p) az
overlay sxxu(1) with sxxuover(ncoef,p)
overlay coeff(1,1) with betadraw(ncoef*p)
overlay vcoefp(1,pk1+1) with gamm0(p,icon)
overlay vcoefp(1,1) with gamma(p,pk1)
overlay vcoefp(1,pk1+1+icon) with alfa(p,2)
overlay eye(1,1) with pull(p,r)
overlay betaptr(1,r+1) with az(p-r,p-r)
overlay mmvec(1,1) with a1(p,p-r)
overlay mmval(1) with d1(p-r,1)
overlay meivec(1,1) with betapcp(p,p-r)
matrix sxx = tr(decomp(xxx))
matrix sinv = tr(decomp(inv(biglm)))
matrix betaols = betadraw
matrix first = second =  const(0.0)
matrix kpswpim = kpswpiv = const(0.0)
scratch(matrix=impblk) nvar*p 1 nstep impulses
list ieqn = 1 to nvar
input betap
 1.0 0.0 -1.0 .1201614
 0.0 1.0 0.0  -1.0
input hld
 0.0 1.0
 1.0 0.0
matrix eye = iden(1.0)
matrix pulh = pull
matrix pltb = tr(pull)*tr(betap)
matrix betac = tr(betap)*inv(pltb)
matrix betacp = tr(betac)
ewise one(i,j) = 1.0
matrix m = eye - tr(betap)*inv(betap*tr(betap))*betap
eigen m meival meivec
overlay meivec(1,1) with betaprp(p,p-r)
do draws = 1,ndraws
cntrl(pause)
display 'iteration = ' draws
matrix ranw = ran(1.0/sqrt(nobs))
matrix rantr = tr(ranw)*sinv
if ndraws>1
     {
     matrix biglm = inv((tr(rantr)*rantr))
     }
matrix swish = decomp(biglm)
matrix ranc = ran(1.0)
matrix sxxuover = tr(sxx)*ranc
if ndraws>1
     {
     matrix betadraw = betaols+kronid(swish,sxxu)
     }
matrix vcoefp = tr(coeff)
matrix gam = gamma
if ndraws==1
     {
     write vcoefp
     display '          Gamma(i) Matricies'
     display ' '
     write gam
     }
matrix alpha = alfa
if ndraws==1
     {
     display '          Alpha Matrix'
     display ' '
     write alpha
     }
matrix bigpi = scale(-1)*alpha*betap
if ndraws==1
     {
     display '           Phi Matrix'
     display ' '
     write bigpi
     }
* Common Trends Analysis and Decompositions
matrix m1rr = inv(biglm)* $
    (eye-alpha*inv(tr(alpha)*inv(biglm)*alpha)*tr(alpha)*inv(biglm))
eigen m1rr mmval mmvec
matrix sd1 = diag(d1)
matrix aprp = a1*decomp(sd1)
if ndraws==1
     {
     display '        Alpha Prp Matrix'
     display ' '
     write aprp
     }
matrix pi11 = eye -gam*kroneker(eye,one) - scale(k-2)*bigpi
matrix c = betaprp*inv(tr(aprp)*pi11*betaprp)*tr(aprp)
if ndraws==1
     {
     display '          C(1) Matrix'
     display ' '
     write c
     }
matrix m = eye - tr(betacp)*inv(betacp*tr(betacp))*betacp
eigen m meival meivec
matrix betaptr = tr(betapcp)
matrix betapcp = tr(inv(az)*betaptr)
if ndraws==1
     {
     display '           Beta Prp Matrix = KPSW A(0)'
     display ' '
     write betapcp
     }
* Matrix to interchange columns of betapcp if necessary
matrix phi = hld
matrix betapcp = betapcp*phi
if ndraws==1
     {
     display '      KPSW A(0) Matrix after Column Exchange'
     display ' '
     write betapcp
     }
matrix f = inv(tr(betapcp)*betapcp)*tr(betapcp)*c
matrix sigeta = f*biglm*tr(f)
if ndraws==1
     {
     display ' Covariance Matrix of Permanent Components'
     display ' '
     write sigeta
     }
matrix phi = decomp(sigeta)
matrix norm1 = iden(1.0)
ieval pr = p-r
do ik = 1,pr
    eval norm1(ik,ik) = 1.0/phi(ik,ik)
    end do ik
matrix phi = phi*norm1
ewise kpswpim(i,j) = kpswpim(i,j) + phi(i,j)
ewise kpswpiv(i,j) = kpswpiv(i,j) + phi(i,j)**2
if ndraws==1
     {
     display '          KPSW PI Matrix'
     display ' '
     write phi
     }
matrix ahat = betapcp*phi
if ndraws==1
     {
     display '     KPSW AHat Matrix'
     display ' '
     write ahat
     }
matrix fstar = inv(tr(ahat)*ahat)*tr(ahat)*c
matrix mm = eye - tr(fstar)*inv(fstar*tr(fstar))*fstar
eigen mm eival g
matrix mm = tr(c)*c
eigen mm eival hinv
matrix h = inv(hinv)
matrix r0 = c*hinv
dimension sigeta(p,p)
matrix r3 = fstar*hinv
matrix r1 = iden(1.0)
do ii = 1, pr
do ij = 1, pr
    eval r1(ii,ij) = r3(ii,ij)
    end do ij
    end do ii
matrix hs = r1*h
matrix hsinv = inv(hs)
matrix sigeta = hs*biglm*tr(hs)
matrix r0 = decomp(sigeta)
matrix r1 = iden(1.0)
do ii = 1, p
    eval r1(ii,ii) = r0(ii,ii)
    end do ii
matrix r0 = inv(r0)
matrix hs = r1*r0*hs
matrix hsinv = inv(hs)
if ndraws==1
     {
     display '      KPSW Gamma(0) Matrix'
     display ' '
     write hsinv
     }
matrix r0 = c*hsinv
if ndraws==1
     {
     display '  Check  C(1)*Gamma(0) =  [Ahat|0]'
     display ' '
     write r0
     }
do ii = 1, p
    impulse(noprint,decomp=hsinv) nvar nstep ii
    cards ieqn impulses+nvar*(ii-1)+ieqn 1
end do ii
ewise first(i,j) = first(i,j)+impblk(i,j)
ewise second(i,j) = second(i,j)+impblk(i,j)**2
matrix sigeta = hs*biglm*tr(hs)
if ndraws== 1
     {
     display ' Covariance Matrix of Eta Disturbances'
     display ' '
     write sigeta
     }
matrix r0 = decomp(sigeta)
matrix hsinvu =hsinv*r0
end do draws
errors(decomp=hsinvu) 12 20
cards ieqn
end
matrix first = scale(1.0/ndraws)*first
matrix second = scale(1.0/ndraws)*second
ewise second(i,j) = second(i,j) - first(i,j)**2
matrix kpswpim = scale(1.0/ndraws)*kpswpim
matrix kpswpiv = scale(1.0/ndraws)*kpswpiv
ewise kpswpiv(i,j) = kpswpiv(i,j) -kpswpim(i,j)**2
open copy a:impulse.out
write(unit=copy) first
write(unit=copy) second
close copy
open copy a:pi.out
write(unit=copy) kpswpim
write(unit=copy) kpswpiv
close copy
end
matrix sigeta = hs*biglm*tr(hs)
display '  Covariance Matrix of Eta Disturbances'
display '       Check: Should be Diagonal'
display ' '
write sigeta
matrix r0 = decomp(sigeta)
matrix hsinvu = hsinv*r0
errors(decomp=hsinvu) 10 40
#1
#2
#3
#4
#5
#6
#7
#8
#9
#10
end

************************************************************************

all 0 40
smpl 1 40
open data a:impulse.unr
declare rect first second
scratch(matrix=first) 48 1 40 ibasf
scratch(matrix=second) 48 1 40 ibass
input(unit=data) first
input(unit=data) second
ieval num1 = 12
ieval num2 = 7
set top1 / = (ibasf+num1)(t) + 1.96*(ibass+num1)(t)
set bot1 / = (ibasf+num1)(t) - 1.96*(ibass+num1)(t)
set top2 / = (ibasf+num2)(t) + 1.96*(ibass+num2)(t)
set bot2 / = (ibasf+num2)(t) - 1.96*(ibass+num2)(t)
grparm(portrait)  subheader 12
spgraph(vfields=2)
graph(samesize,header='Figure A1',number=0, $
   subhead='Velocity irf to Money Growth Shock') 3
#ibasf+num1  / 1
#top1 / 2
#bot1 /2
print / bot1 ibasf+num1 top1
graph(samesize,header='Figure A2',number=0, $
   subhead='Real Output irf to Money Growth Shock') 3
#ibasf+num2  / 1
#top2 / 2
#bot2 /2
print / bot2 ibasf+num2 top2
spgraph(done)
ieval num1 = 11
ieval num2 = 6
set top1 / = (ibasf+num1)(t) + 1.96*(ibass+num1)(t)
set bot1 / = (ibasf+num1)(t) - 1.96*(ibass+num1)(t)
set top2 / = (ibasf+num2)(t) + 1.96*(ibass+num2)(t)
set bot2 / = (ibasf+num2)(t) - 1.96*(ibass+num2)(t)
spgraph(vfields=2)
graph(samesize,header='Figure A3',number=0, $
   subhead='Inflation irf to Money Growth Shock') 3
#ibasf+num1  / 1
#top1 / 2
#bot1 /2
print / bot1 ibasf+num1 top1
graph(samesize,header='Figure A4',number=0, $
   subhead='Money Growth irf to Money Growth Shock') 3
#ibasf+num2  / 1
#top2 / 2
#bot2 /2
print / bot2 ibasf+num2 top2
spgraph(done)
ieval num1 = 8
ieval num2 = 24
set top1 / = (ibasf+num1)(t) + 1.96*(ibass+num1)(t)
set bot1 / = (ibasf+num1)(t) - 1.96*(ibass+num1)(t)
set top2 / = (ibasf+num2)(t) + 1.96*(ibass+num2)(t)
set bot2 / = (ibasf+num2)(t) - 1.96*(ibass+num2)(t)
spgraph(vfields=2)
graph(samesize,header='Figure A5',number=0, $
   subhead='T-bill Rate irf to Money Growth Shock') 3
#ibasf+num1  / 1
#top1 / 2
#bot1 /2
print / bot1 ibasf+num1 top1
graph(samesize,header='Figure A6',number=0, $
   subhead='Velocity irf to Real Output Shock') 3
#ibasf+num2  / 1
#top2 / 2
#bot2 /2
print / bot2 ibasf+num2 top2
spgraph(done)
ieval num1 = 19
ieval num2 = 23
set top1 / = (ibasf+num1)(t) + 1.96*(ibass+num1)(t)
set bot1 / = (ibasf+num1)(t) - 1.96*(ibass+num1)(t)
set top2 / = (ibasf+num2)(t) + 1.96*(ibass+num2)(t)
set bot2 / = (ibasf+num2)(t) - 1.96*(ibass+num2)(t)
spgraph(vfields=2)
graph(samesize,header='Figure A7',number=0, $
   subhead='Real Output irf to Real Output Shock') 3
#ibasf+num1  / 1
#top1 / 2
#bot1 /2
print / bot1 ibasf+num1 top1
graph(samesize,header='Figure A8',number=0, $
   subhead='Inflation irf to Real Output Shock') 3
#ibasf+num2  / 1
#top2 / 2
#bot2 /2
print / bot2 ibasf+num2 top2
spgraph(done)
ieval num1 = 18
ieval num2 = 20
set top1 / = (ibasf+num1)(t) + 1.96*(ibass+num1)(t)
set bot1 / = (ibasf+num1)(t) - 1.96*(ibass+num1)(t)
set top2 / = (ibasf+num2)(t) + 1.96*(ibass+num2)(t)
set bot2 / = (ibasf+num2)(t) - 1.96*(ibass+num2)(t)
spgraph(vfields=2)
graph(samesize,header='Figure A9',number=0, $
   subhead='Money Growth irf to Real Output Shock') 3
#ibasf+num1  / 1
#top1 / 2
#bot1 /2
print / bot1 ibasf+num1 top1
graph(samesize,header='Figure A10',number=0, $
   subhead='T-bill Rate irf to Real Output Shock') 3
#ibasf+num2  / 1
#top2 / 2
#bot2 /2
print / bot2 ibasf+num2 top2
spgraph(done)
end