OUTPUT FILE = RESULT.OUT RESET;
FORMAT /m1/ro  12,4;
@--------------------------------------------------------------------@
@ Set Autoregressive Parameters                                      @
@--------------------------------------------------------------------@
OUTPUT OFF;
nvar = 2;

nlag = 2;

Let auto[2,4]      =   -.20   .00     -.01   .00
                        .00  -.20      .00  -.01  ;

Let sigma[2,2] =        .01  .00  .00   .01 ;

Let b[2,1]  =    0 0 ;

Let nval[1,2] =   3  3 ;
@--------------------------------------------------------------------@
@ Set parameters for asset pricing                                   @
@--------------------------------------------------------------------@
beta   =  .97 ;   @ subjective discount factor                       @
psi    = 1.30 ;   @ coefficient of relative risk aversion            @
convbl = 1;       @ variable number that is log consumption growth   @
assets = 2;       @ vector of variable numbers of assets to be       @
                  @     priced                                       @
hold   = 1;       @ holding period for assets                        @
maxit  = 1000;    @ maximum iterations in computing asset prices     @
nforward = 4;     @ number of periods forward for which prices of    @
                  @    pure discount bills are to be computed        @
@--------------------------------------------------------------------@
@ Load data array                                                    @
@--------------------------------------------------------------------@
Load  ghmat[55,3] = ghquad.dat;
@--------------------------------------------------------------------@
@ Clear procs and variables with global references                   @
@--------------------------------------------------------------------@
Proc pisub; retp(""); Endp;
Proc pifn; retp(""); Endp;
Proc reachrow; retp(""); Endp;
Proc hzdirpr; retp(""); Endp;
Clear ns, nstate, cntparms, exval, pimat, pista, ymat, wvec, reach;
@--------------------------------------------------------------------@
@ Define additional parameters and print parameters                  @
@--------------------------------------------------------------------@
arparms = auto~b~sigma;
siginv = inv(sigma);
detsigma = det(sigma);
ns = prodc(nval');
nstate = ns^nlag;
OUTPUT ON;
print "Autoregressive Parameters";
dummy = prmat(auto,"auto");
dummy = prmat(b,"b");
dummy = prmat(sigma,"sigma");
print;
print "Control Parameters"; print;
print /m1/rz "     nval   = " nval;
print /m1/rz "     nlag   = " nlag;
print /m1/rz "     ns     = " ns;
print /m1/rz "     nstate = " nstate;
print;
print "Asset Pricing Parameters"; print;
print "     beta   = " beta;
print "     psi    = " psi ;
OUTPUT OFF;
@--------------------------------------------------------------------@
@ Calculate variance/covariance matrix of stationary distribution    @
@--------------------------------------------------------------------@
sigsta = arsta(auto,sigma,1);
OUTPUT ON;
dummy = prmat(sigsta[1:nvar,1:nvar],"sigsta");
print;
OUTPUT OFF;
@--------------------------------------------------------------------@
@ Calculate markov chain                                             @
@--------------------------------------------------------------------@
dummy = mchain(auto,b,sigma,nlag,nval,ghmat,&armean,&code,&decode,
               &hzdirpr,&markstat,&pisub,&pifn,&reachrow,&setquad,
               &z_to_y) ;
@--------------------------------------------------------------------@
@ Compare induced autoregressive paramaters                          @
@--------------------------------------------------------------------@
print;
print "Calculating the linear projection of the variables on their lags using
the discrete laws of motion";
indparms = arproj(pimat,pista,ymat,&decode,&hzdirpr);
OUTPUT ON;
print;
print "Induced Autoregressive Parameters";
dummy = prmat(indparms[.,1:nvar*nlag],"auto_ind");
dummy = prmat(indparms[.,nvar*nlag+1],"b_ind");
dummy = prmat(indparms[.,nvar*nlag+2:cols(indparms)],"sig_ind");
OUTPUT OFF;
@--------------------------------------------------------------------@
@ Set up marginal rate of substitution matrix;                       @
@--------------------------------------------------------------------@
congr = ymat[.,convbl];
mrsmat = crr_mrs( ones(nstate/ns,1).*.ymat );
@--------------------------------------------------------------------@
@ Price the assets.  The array v contains the price/dividend ratios  @
@--------------------------------------------------------------------@
divgr = exp( ymat[.,assets] ) ~ ones(rows(ymat),1);
prntrace = 1;
v = prasset0(maxit,divgr,prntrace,&mexpand,);
OUTPUT ON;
print;
print "The price/dividend ratios v are";
print v;
OUTPUT OFF;
@--------------------------------------------------------------------@
@ derive the term structure of interest rates                        @
@--------------------------------------------------------------------@
ftilde = ones(nstate,ns);
i_forw = 1;
term = zeros(nstate,nforward);
Do while (i_forw le nforward);
   term[.,i_forw] = sumc((pimat.*mrsmat.*ftilde)');
   ftilde = reshape(term[reach,i_forw],nstate,ns);
   i_forw = i_forw + 1;
Endo;
OUTPUT ON;
print;
print "The prices of pure discount bills are";
print term;
print;
OUTPUT OFF;
#INCLUDE "gprocs.gpr"
END;