@--------------------------------------------------------------------@ @ Procedures @ @--------------------------------------------------------------------@ @--------------------------------------------------------------------@ @ code to load exec's @ @--------------------------------------------------------------------@ hdpc=zeros(91,1); @ allocate buffer for code @ loadexe hdpc=hdp.gxe; @ load function into buffer @ @.hr proc code @ @@ @--------------------------------------------------------------------@ @ proc code @ @--------------------------------------------------------------------@ Proc Code(imat,base); Local nr, nc; @--------------------------------------------------------------------@ @ Last Modified: March 31, 1986 10:00 am @ @ Written: March 15, 1986 8:15 pm @ @--------------------------------------------------------------------@ @--------------------------------------------------------------------@ @ Proc to code each row of imat as @ @ @ @ x = n1 + (n2-1)*base + ... + (nL-1)*base(L-1) @ @ @ @ where (n1,n2,...,nL) is a row of imat and L = cols(imat). @ @ @ @ Input: imat matrix @ @ base scaler @ @ Returns: vector rows(imat)x1 @ @--------------------------------------------------------------------@ @--------------------------------------------------------------------@ @ Check for errors @ @--------------------------------------------------------------------@ If ( cols(base) gt 1 or rows(base) gt 1 or base[1,1] lt 1 ); Print " *** Code error: BAD BASE " ; Print " base = "; Print base; Print " imat = "; Print imat; retp(""); Endif; If ( maxc( maxc(imat) ) gt base or maxc( maxc(-imat)) ge 0 ); Print " *** Code error: An element of the input array either exceeds "; Print " base or is less than or equal to zero"; Print; Print " input array = "; Print imat; Print " base = "; Print base; retp(""); Endif; @--------------------------------------------------------------------@ @ Start main part of proc @ @--------------------------------------------------------------------@ nr = rows(imat); nc = cols(imat); If (nc eq 1); retp(imat); Else; retp( ( imat - (zeros(nr,1)~ones(nr,nc-1)) )*seqm(1,base,nc) ); Endif; Endp; @ Save code; @ @.hr proc decode @ @@ Proc decode(x,base,L); Local nr; @--------------------------------------------------------------------@ @ Last Modified: March 17, 1986 10:15 am @ @ Written: March 15, 1986 11:55 am @ @--------------------------------------------------------------------@ @--------------------------------------------------------------------@ @ Proc to decode the rows of x where a typical row is of the form @ @ @ @ x[k,.] = n1 + (n2-1)*base^1 + ... (nL-1)*base^(L-1) @ @ @ @ Input: x vector rows(x)x1 @ @ L scaler @ @ Returns: matrix rows(x)xL @ @--------------------------------------------------------------------@ @--------------------------------------------------------------------@ @ Check for errors @ @--------------------------------------------------------------------@ If ( base le 0 or L lt 1); Print "*** Decode error bad base = " base " or bad length = " L; retp(""); Endif; If (cols(x) gt 1); Print "*** Decode error: input array is not a column vector"; Print " input array = "; Print x; retp(""); Endif; If ( maxc( maxc(-x) ) ge 0 ); Print "*** Decode warning: at least one element of input vector"; Print " is less than or equal to zero "; Print " input array = "; Print x; Endif; @--------------------------------------------------------------------@ @ Start proc @ @--------------------------------------------------------------------@ nr = rows(x); x = x - ones(nr,1); retp( (trunc( x./(seqm(1,base,L))' )) % base + ones(nr,L) ); Endp; @ Save decode; @ @.hr proc mexpand @ @@ @--------------------------------------------------------------------@ @ proc mexpand @ @--------------------------------------------------------------------@ proc mexpand(a,index); Local indexmin, indexmax, nc, b, m, k, bcol, indexcol ; @--------------------------------------------------------------------@ @ procedure to "expand" the matrix a according to the matrix index @ @ @ @ returns: a[index[i,j],m] @ @--------------------------------------------------------------------@ @--------------------------------------------------------------------@ @ Initialize and check for errors @ @--------------------------------------------------------------------@ indexmax = maxc(maxc(index)); indexmin = -maxc(maxc(-index)); If ( indexmin lt 1 or indexmax gt rows(a) ); Print; Print /m1/rz "*** mexpand error: an index is out of range"; Print; Print /m1/rz "argument 1 is " rows(a) "x" cols(a); Print /m1/rz "argument 2 is " rows(index) "x" cols(index); Print; Print /m1/rz " smallest element of argument 2(index) is " indexmin; Print /m1/rz " largest element of argument 2(index) is " indexmax; Print /m1/rz " rows of argument 1 = " rows(a) ; Print; retp(""); Endif; nc = cols(index); b = zeros(rows(index),cols(a)*cols(index)); @--------------------------------------------------------------------@ @ Start main part of proc @ @--------------------------------------------------------------------@ m = 1; Do while ( m le cols(a) ); k = 1; Do while (k le nc); indexcol = index[.,k]; bcol = (m-1)*nc + k; b[.,bcol] = a[indexcol,m]; k = k + 1; Endo; m = m + 1; Endo; retp(b); Endp; @.hr proc phimat @ @@ Proc phimat(xmat,mu,siginv,detsigma); Local zmat; @--------------------------------------------------------------------@ @ Proc for row-wise multvariate normal distribution @ @--------------------------------------------------------------------@ @--------------------------------------------------------------------@ @ Input: xmat matrix n x nvar @ @ mu vector 1 x nvar @ @ siginv matrix nvar x nvar @ @ detsigma scaler @ @ @ @ Returns: vector n x 1 @ @ @ @ Notes: The elements of returned vector are the multivariate @ @ normal density N(mu,sigma) evaluated at each row of xmat @ @ @ @--------------------------------------------------------------------@ zmat = sumc( ( ((xmat - mu)*siginv).*(xmat - mu) )' ); retp( exp( -0.5*zmat )/sqrt(2*pi*detsigma) ); Endp; @ Save phimat; @ @.hr proc prmat @ @@ @--------------------------------------------------------------------@ @ Proc to print a matrix @ @--------------------------------------------------------------------@ proc prmat(x,name); Print; If ( rows(x) gt 4 or cols(x) gt 4 ); Print /m1/rz name "[" rows(x) "," cols(x) "] = "; Else; Print /m1/rz name " ="; Endif; Print /m1/ro x; retp(""); Endp; @ Save prmat; @ @.hr proc armean @ @@ Proc armean(auto,b); Local nvar, nlag, temp, k; @--------------------------------------------------------------------@ @ Last modified: March 19, 1986 11:15 am @ @ Written: March 19, 1986 11:15 am @ @--------------------------------------------------------------------@ @ Proc to calculate the mean of an autoregressive process @ @ @ @ Input: auto matrix nrxnc @ @ b vector nrx1 @ @ Returns mean vector nrx1 @ @ @ @--------------------------------------------------------------------@ @--------------------------------------------------------------------@ @ Check for errors @ @--------------------------------------------------------------------@ If ( rows(auto) ne rows(b) or trunc( cols(auto)/rows(auto) ) ne cols(auto)/rows(auto) ); Print " *** armean error: nonconformable input data "; Print on; auto=auto; b=b; Print off; retp(""); Endif; @--------------------------------------------------------------------@ @ Start main part of proc @ @--------------------------------------------------------------------@ nvar = rows(auto); nlag = cols(auto)/rows(auto); temp = zeros(nvar,nvar); k = 1; Do while (k le nlag); temp = temp + auto[.,( 1+ (k-1)*nvar ):k*nvar]; k = k + 1; Endo; retp( (Inv( eye(nvar) - temp )*b) ); Endp; @ save armean; @ @.hr proc markstat @ @@ Proc markstat(pimat,code,decode); Local pista, temp, epsta, maxstit, crit, numstit, j, k, st, stold, nstate, ns, nlag, base, jold, delta, code:proc, decode:proc ; @--------------------------------------------------------------------@ @ Last modified: March 19, 1986 11:15 am @ @ Written: March 19, 1986 11:15 am @ @--------------------------------------------------------------------@ @ Proc to calculate the stationary distribution of the markov chain @ @ @ @ Input: pimat matrix nstate x ns @ @ code proc @ @ decode proc @ @ Returns: pista vector nstate x 1 @ @ @ @--------------------------------------------------------------------@ @--------------------------------------------------------------------@ @ Set array sizes and check for errors @ @--------------------------------------------------------------------@ nstate = rows(pimat); ns = cols(pimat); nlag = ln(nstate)/ln(ns); base = ns; @--------------------------------------------------------------------@ @ intitialize @ @--------------------------------------------------------------------@ pista = ones(nstate,1)/nstate; temp = zeros(nstate,1); epsta = .0001; maxstit = 50; crit = 10000; print; numstit = 1; @--------------------------------------------------------------------@ @ Main loop to compute the chain's stationary distribution @ @--------------------------------------------------------------------@ Do while (crit gt epsta and numstit le maxstit); If (nlag eq 1); temp = (pista'*pimat)'; Else; j=1; Do while ( j le nstate ); st = decode(j,base,nlag); stold = ( st[1,2:nlag].*ones(ns,1) ) ~ seqa(1,1,ns) ; jold = code(stold,base); temp[j,1] = pimat[jold,st[1,1]]'*pista[jold,1]; j=j+1; Endo; Endif; delta = abs(temp - pista); crit = maxc(delta); pista = temp; Print " numstit = " numstit "crit = " crit; numstit = numstit + 1; Endo; If ( numstit gt maxstit ); Print "*** Unsuccessful iterations for stationary distribution"; Print on; numstit pista delta; Endif; retp(pista); Endp; @ Save markstat; @ @.hr proc setquad @ @@ @--------------------------------------------------------------------@ @ proc setquad @ @--------------------------------------------------------------------@ Proc setquad(nval,nvar,quadat); Local np, qmat, zmat, wvec, m, k, k1 ; @--------------------------------------------------------------------@ @ Procedure to put all possible combinations of quadrature abscissa @ @ and corresponding products of quadrature weights in @ @ returned matrix @ @--------------------------------------------------------------------@ @ Last Modified: March 22, 1986 10:40 am @ @ Written: March 22, 1986 10:40 am @ @ @ @--------------------------------------------------------------------@ @--------------------------------------------------------------------@ @ Start main procedure @ @--------------------------------------------------------------------@ np = nval[1,1]; Gosub SELECTQD; zmat = qmat[.,1]; wvec = qmat[.,2]; If (nvar ge 2 ); m = 2; Do while (m le nvar); np = nval[1,m]; Gosub SELECTQD; zmat = ( ones(np,1).*.zmat ) ~ ( qmat[.,1].*.ones(rows(zmat),1) ); wvec = qmat[.,2].*.wvec; m = m + 1; Endo; Endif; retp(zmat~wvec); @--------------------------------------------------------------------@ @ Subroutine to select from quadat the absciss and weights @ @ for np point quadrature rule @ @--------------------------------------------------------------------@ SELECTQD: k1 = ( (np-1).*np )/2 + 1; k = k1; qmat = zeros(np,2); do while ( k le k1-1 + np ); qmat[k-k1+1,1:2] = quadat[k,2:3]; k=k+1; endo; Return; Endp; @ Save setquad; @ @.hr proc z_to_y @ @@ proc z_to_y(zmat,exval,sigma); @--------------------------------------------------------------------@ @ Last modified: March 19, 1986 11:15 am @ @ Written: March 19, 1986 11:15 am @ @--------------------------------------------------------------------@ @ Proc to transform the z's to the y's @ @ @ @ Input: ymat matrix nsxnvar @ @ sigma matrix nvarxnvar @ @ exval vector 1xnvar @ @ Returns: matrix nsxnvar @ @ @ @--------------------------------------------------------------------@ @--------------------------------------------------------------------@ @ Check for errors @ @--------------------------------------------------------------------@ If (cols(zmat) ne cols(sigma) or cols(sigma) ne rows(sigma) or cols(exval) ne cols(zmat) ); Print "*** z_to_y error: nonconformable input matrices"; Print on; zmat=zmat; sigma=sigma; Print off; retp(""); Endif; @--------------------------------------------------------------------@ @ Start main part of proc @ @--------------------------------------------------------------------@ retp( zmat*chol(sigma) + ones(rows(zmat),1).*.exval ) ; Endp; @.hr proc mchain @ @@ @--------------------------------------------------------------------@ @ proc mchain @ @--------------------------------------------------------------------@ proc mchain(auto,b,sigma,nlag,nval,ghmat, armean, code, decode, hzdirpr, markstat, pisub, pifn, reachrow, setquad, z_to_y) ; Local armean:proc, code:proc, decode:proc, hzdirpr:proc, markstat:proc, pisub:proc, pifn:proc, reachrow:proc, setquad:proc, z_to_y:proc; Local nvar, zwmat, zmat, wvec, ns, nstate, cntparms, base, j, is, reach, err, exval, ymat, pimat, pista; @--------------------------------------------------------------------@ @ procedure to calculate discrete markov chain that approximates VAR @ @--------------------------------------------------------------------@ @--------------------------------------------------------------------@ @ Initialize and check for errors @ @--------------------------------------------------------------------@ nvar = rows(auto); If (cols(auto) ne nvar*nlag); Print "Auto assignment error"; Print nvar nlag auto; Retp(""); Endif; If (rows(sigma) ne nvar or cols(sigma) ne nvar); Print "Sigma assigment error"; Print nvar sigma; Retp(""); Endif; If (rows(b) ne nvar); Print "b assigment error"; Print nvar b; Retp(""); Endif; @--------------------------------------------------------------------@ @ Start of Main part of proc @ @--------------------------------------------------------------------@ @--------------------------------------------------------------------@ @ Set quadrature arrays @ @--------------------------------------------------------------------@ zwmat = setquad(nval,nvar,ghmat); zmat = zwmat[.,1:(cols(zwmat)-1)]; wvec = zwmat[.,cols(zwmat)]; err = Varput(wvec,"wvec"); If (err eq 0); gosub vputerr; Endif; Clear zwmat; @--------------------------------------------------------------------@ @ set additional constants and arrays of parameters @ @--------------------------------------------------------------------@ ns = rows(zmat); nstate = ns^nlag; cntparms = nvar|nlag|ns|nstate; base = ns; err = Varput(ns,"ns"); If (err eq 0); gosub vputerr; Endif; err = Varput(nstate,"nstate"); If (err eq 0); gosub vputerr; Endif; err = Varput(cntparms,"cntparms"); If (err eq 0); gosub vputerr; Endif; @--------------------------------------------------------------------@ @ Form the reachable states and put them in reach @ @--------------------------------------------------------------------@ Print "Forming the reachable states"; reach = ones(nstate,ns); j = 1; Do while ( j le nstate ); is = decode(j,base,nlag); reach[j,.] = reachrow(is,ns,base); j = j + 1; Endo; print; print "Reachable states formed"; err = Varput(reach,"reach"); If (err eq 0); gosub vputerr; Endif; @--------------------------------------------------------------------@ @ Calculate the mean of the process @ @--------------------------------------------------------------------@ exval = armean(auto,b)'; err = Varput(exval,"exval"); If (err eq 0); gosub vputerr; Endif; @--------------------------------------------------------------------@ @ Transform to the Y's and add in the means @ @--------------------------------------------------------------------@ Print on; ymat = z_to_y(zmat,exval,sigma); clear zmat; Print off; err = Varput(ymat,"ymat"); If (err eq 0); gosub vputerr; Endif; @--------------------------------------------------------------------@ @ Get transition probabilities @ @--------------------------------------------------------------------@ Print; Print "Starting to compute transition probabilities"; pimat = pisub( seqa(1,1,nstate), ymat, ns, ymat, ns ); Print; Print "Transition probabilities computed"; err = Varput(pimat,"pimat"); If (err eq 0); gosub vputerr; Endif; @--------------------------------------------------------------------@ @ Compute the chain's stationary distribution @ @--------------------------------------------------------------------@ pista = markstat(pimat,&code,&decode); Print; Print "Stationary distribution computed"; err = Varput(pista,"pista"); Retp(""); vputerr: Print; Print "**** Mchain error -- symbol table full. Variables not set"; END; Return; Endp; @ Save mchain; @ @.hr proc prassset @ @@ @--------------------------------------------------------------------@ @ proc prasset @ @--------------------------------------------------------------------@ @ dummy defines of pimat, mrsmat, and reach @ pimat = 1; mrsmat = 1; reach = 1; proc prasset0(maxit,divgr,prntrace,mexpand); Local mexpand:proc, nasset, nstate, v, numit, mincrit, m, crit, mrsbig, pibig, i, pitsum, pitilde, temp, nextv ; @--------------------------------------------------------------------@ @ procedure to calculate equilibrium asset prices @ @--------------------------------------------------------------------@ @--------------------------------------------------------------------@ @ intialize @ @--------------------------------------------------------------------@ nasset = cols(divgr); nstate = rows(pimat); v = zeros(nstate,nasset); numit = ones(nasset,1); mincrit = .0001; Print; Print "Pricing the assets"; If ( ( rows(reach) ne rows(pimat) ) or ( cols(reach) ne cols(pimat) ) ); Print " #### Warning -- reach has been reset. If prassets iterates "; Print " it will crash. "; Print; Endif; @--------------------------------------------------------------------@ @ start main part of procedure @ @--------------------------------------------------------------------@ save reach; reach = 0; m = 1; Do while ( m le nasset ); pitilde = pimat.*mrsmat.*divgr[.,m]'; save pimat; save mrsmat; pimat = 1; mrsmat = 1; pitsum = sumc(pitilde'); load reach; If ( nstate^2 le 8100 ); pibig = zeros(nstate,nstate); i = 1; do while (i le nstate); pibig[i,reach[i,.]] = -pitilde[i,.]; pibig[i,.] = (seqa(1,1,nstate)-i .== 0)' + pibig[i,.]; i = i + 1; endo; pitilde = 0; reach = 0; v[.,m] = (Inv( pibig ))*pitsum; Else; @ Iterate @ crit = mincrit + 1; Do while ( crit ge mincrit and numit[m,1] le maxit ); nextv = sumc( ( pitilde.*(1 + mexpand(v[.,m],reach)) )' ); crit = maxc( abs(v[.,m] - nextv) ); v[.,m] = nextv; If (prntrace eq 1 and numit[m,1]%10 eq 1); Print; Print /m1/rz "asset = " m "crit = " crit "numit = " numit[m,1]; If (numit[m,1]%50 eq 1 or crit le 10*mincrit ); Print; Print /m1/rz nextv'; Endif; Endif; numit[m,1] = numit[m,1] + 1; Endo; pitilde = 0; reach = 0; Endif; load pimat; load mrsmat; m = m + 1; Endo; load reach; If ( maxc(numit) gt maxit and maxit gt 1) ; Print; Print " Asset pricing iterations failed to converge for at least "; Print " one of the assets. "; Print; Print /m1/rz "maximum iterations allowed = " maxit; Print; Print /m1/rz "iterations attempted for each asset are : "; Print; Print; /m1/rz (seqa(1,1,nasset))~numit; Print; Endif; Print; Print "Assets priced"; retp(v); Endp; @.hr proc arsta @ @@ @--------------------------------------------------------------------@ @ proc arsta @ @--------------------------------------------------------------------@ proc arsta(auto,sigma,prntrace); Local nvar, nlag, eps, maxit, mincrit, nn, arbig, crit, numit, temp, sigsta; @--------------------------------------------------------------------@ @ Procedure to calculate the autoregressions stationary distribution @ @--------------------------------------------------------------------@ @--------------------------------------------------------------------@ @ Check for errors @ @--------------------------------------------------------------------@ nvar = rows(auto); nlag = cols(auto)/nvar; If ( abs( round(nlag) - nlag ) ge 1.0e-5 or rows(sigma) ne nvar or cols(sigma) ne nvar ); Print; Print; " *** arsta error: nonconformable input data ";Print; Print " auto is " rows(auto) "x" cols(auto); Print " sigma is " rows(sigma) "x" cols(sigma); Print; Print " auto = "; Print auto; Print; Print " sigma ="; Print sigma; Print; retp(""); Endif; @--------------------------------------------------------------------@ @ Initialize @ @--------------------------------------------------------------------@ eps = .0001; maxit = 200; mincrit = eps*maxc(maxc( abs(sigma) )); @--------------------------------------------------------------------@ @ start main part of proc; @ @--------------------------------------------------------------------@ nn = nlag-1; If ( nn eq 0); arbig = auto; sigsta = sigma; Else; arbig = auto | ( eye(nn*nvar)~zeros(nn*nvar,nvar) ); sigsta = diagrv( zeros(nlag,nlag),1|zeros(nn,1) ).*.sigma; Endif; print; print "Calculating stationary error covariance matrix"; print; numit = 1; crit = mincrit + 1; temp = sigsta; Do while ( crit ge mincrit and numit le maxit); temp = arbig*temp*arbig'; sigsta = sigsta + temp; crit = maxc(maxc( abs(temp) )); If (prntrace eq 1); Print " Iteration for sigsta: numit = " numit "crit = " crit; Endif; numit = numit + 1; Endo; If (numit le maxit); retp(sigsta); @ done @ Else; @ otherwise unsuccesful @ Print; Print "Unsuccesful iterations for stationary covariance matrix "; Print; Print "numit = " numit "maxit = " maxit; Print "crit = " crit "mincrit = " mincrit; Print; Print "auto = "; Print auto; Print; Print "sigma = "; Print sigma; Print; retp(""); Endif; Endp; @ Save arsta; @ @.hr proc crr_mrs @ @@ proc crr_mrs(xmat); Local ymat, psi, beta, congr, mrsmat; ymat = varget("ymat"); congr = varget("congr"); beta = varget("beta"); psi = varget("psi"); mrsmat = beta*( exp(-psi*congr') ).*ones(rows(xmat),1); retp(mrsmat); Endp; @.hr proc arproj @ @@ proc arproj(pimat,pista,ymat,decode,hzdirpr); Local maxnum, jstep, kmat, yhat, m; Local ns, nstate, nvar, nlag, ybar, ex_yy, ex_yylag; Local j1, jvec, ylag; Local seqnv, syy0, syy, j, k, idj, idk, idkj; Local autoproj, apsum, L, b0, sigbar; Local decode:proc, hzdirpr:proc; @--------------------------------------------------------------------@ @ proc to calculate the linear projection of y on lagged y's using @ @ the discrete laws of motion @ @--------------------------------------------------------------------@ @--------------------------------------------------------------------@ @ initialize and check for errors @ @--------------------------------------------------------------------@ ns = cols(pimat); nstate = rows(pimat); nvar = cols(ymat); If ( (nstate gt ns) and (ns gt 1) ); nlag = ln(nstate)/ln(ns); Else; nlag = 1; Endif; If ( nstate lt ns ); Print; Print " *** arproj error: rows(pimat) < cols(pimat) "; Print; Print " *** exiting procedure "; retp(""); Endif; If ( rows(ymat) ne ns ); Print; Print " *** arproj error: rows(ymat) /= cols(pimat) "; Print; Print " *** exiting procedure "; retp(""); Endif; @--------------------------------------------------------------------@ @ calculate jstep so as to keep the size of matrices within bounds @ @--------------------------------------------------------------------@ maxnum = 8100; jstep = nstate; Do While ( jstep*(nvar^2)*nlag ge maxnum ); jstep = jstep/ns; Endo; @--------------------------------------------------------------------@ @ form conditional expectation @ @--------------------------------------------------------------------@ Print; Print " Forming conditional mean of discrete process"; kmat = ones(nstate,1).*seqa(1,1,ns)'; yhat = zeros(nstate,nvar); m = 1; Do While ( m le nvar ); yhat[.,m] = sumc((pimat.* reshape(ymat[kmat,m],nstate,ns))'); m = m + 1; Endo; ybar = pista'*yhat; @--------------------------------------------------------------------@ @ compute expected values of cross products @ @--------------------------------------------------------------------@ Print; Print " Forming expected values of cross products" ; ex_yylag = zeros(1,(nvar^2)*nlag); ex_yy = zeros(1,nvar^2); j1 = 1; Do while ( j1 le nstate ); jvec = seqa(j1,1,jstep); kmat = decode(jvec,ns,nlag); ylag = reshape(ymat[kmat,.],rows(kmat),nvar*nlag); ex_yy = ex_yy + pista[jvec,.]'* hzdirpr(ymat[kmat[.,1],.]-ybar,ymat[kmat[.,1],.]); ex_yylag = ex_yylag + pista[jvec,.]'* hzdirpr(yhat[jvec,.]-ybar,ylag); j1 = j1 + jstep; Endo; ex_yy = reshape(ex_yy,nvar,nvar); ex_yylag = reshape(ex_yylag,nvar,nvar*nlag); @--------------------------------------------------------------------@ @ form time series matrix of all expected cross products @ @--------------------------------------------------------------------@ seqnv = seqa(1,1,nvar); syy0 = ex_yy~ex_yylag; syy = zeros(nvar*(nlag+1),nvar*(nlag+1)); j = 0; Do while ( j le nlag ); k = j; Do while ( k le nlag ); idj = j*nvar+seqnv; idk = k*nvar+seqnv; idkj = (k-j)*nvar+seqnv; syy[idj,idk] = syy0[seqnv,idkj]; If ( k gt j ); syy[idk,idj] = syy0[seqnv,idkj]'; Endif; k = k +1; Endo; J = j + 1; Endo; @--------------------------------------------------------------------@ @ compute coefficients of linear projection @ @--------------------------------------------------------------------@ Print; Print " Calculating projection parameters"; autoproj = syy[1:nvar,(1+nvar):cols(syy)]* inv( syy[(1+nvar):cols(syy),(1+nvar):cols(syy)] ); apsum = zeros(nvar,nvar); L = 1; Do while ( L le nlag ); apsum = apsum + autoproj[.,(L-1)*nvar+1:L*nvar]; L = L + 1; Endo; @--------------------------------------------------------------------@ @ compute intercepts @ @--------------------------------------------------------------------@ b0 = (eye(nvar) - apsum)*ybar'; @--------------------------------------------------------------------@ @ compute residual covariance matrix @ @--------------------------------------------------------------------@ sigbar = syy[1:nvar,1:nvar] - autoproj* syy[nvar+1:rows(syy),nvar+1:rows(syy)]*autoproj'; @--------------------------------------------------------------------@ @ return results @ @--------------------------------------------------------------------@ retp(autoproj~b0~sigbar); Endp;