*
* Observable index model
* From Sargent & Sims(1977), "Business cycle modeling without pretending
* to have too much a priori economic theory," in "New Methods in
* Business Cycle Research: Proceedings from a Conference", Federal
* Reserve Bank of Minnepolis.
*
* This differs somewhat from the model actually used in the paper above
* because we allow for contemporaneous loadings on the indexes. The model is
*
*    Y(t)=A(L)C(L)Y(t)+D(L)Y(t-1)+v(t)
*
* where Z(t)=C(L)Y(t) are the "indexes", D(L) is diagonal and v(t) is
* assumed to be uncorrelated across variables. Thus, all interactions
* among the variables are captured in the A(L)Z(t) term. The model takes
* the form of a restricted VAR. With dim(Z)<<dim(Y), it will generally
* have many fewer free parameters than the corresponding full VAR.
*
* The A(L) and C(L) aren't separately identifiable without restrictions.
* They obviously need to be normalized for scale and rotation since
* A(L)C(L)=A(L)GG^-1C(L). It's also possible to have a dynamic
* invertible factor that cancels, that is, A(L)C(L)=A(L)G(L)G(L)^-1C(L),
* so a normalization which fixes the lags of C is also necessary.
*
open data bbedata.rat
calendar(m) 1959
data(format=rats) 1959:01 2001:07 ipp ipf ipc ipcd ipcn ipe ipi ipm ipmd ipmnd ipmfg ipd ipn ipmin iput $
 ip ipxmca pmi pmp gmpyq gmyxpq lhel lhelx lhem lhnag lhur lhu680 lhu5 lhu14 lhu15 lhu26 lpnag lp lpgd $
 lpmi lpcc lpem lped lpen lpsp lptu lpt lpfr lps lpgov lphrm lpmosa pmemp gmcq gmcdq gmcnq gmcsq gmcanq $
 hsfr hsne hsmw hssou hswst hsbr hmob pmnv pmno pmdel mocmq msondq fsncom fspcom fspin fspcap fsput fsdxp $
 fspxe exrsw exrjan exruk exrcan fyff fygm3 fygm6 fygt1 fygt5 fygt10 fyaaac fybaac sfygm3 sfygm6 sfygt1 $
 sfygt5 sfygt10 sfyaaac sfybaac fm1 fm2 fm3 fm2dq fmfba fmrra fmrnba fclnq fclbmc ccinrv pmcp pwfsa pwfcsa $
 pwimsa pwcmsa psm99q punew pu83 pu84 pu85 puc pucd pus puxf puxhs puxm lehcc lehm hhsntn
*
* NLAGSC is the number of lags in C. It includes the 0 lag as well.
* NLAGSD is the number of lags in D.
* NLAGSA is the number of lags in A. It includes the 0 lag as well.
* NLAGSRF is the derived number of lags in the reduced form.
*
compute nlagsc=2
compute nlagsd=2
compute nlagsa=2
compute nlagsrf=%imax(nlagsc+nlagsa,nlagsd)
*
* NF is the number of factors (indexes)
*
compute nf=1
*
* Specifies the reduced form for the VAR. The variable list is taken
* from this.
*
system(model=reducedform)
variables ip ipd ipn ipp ipc ipcd ipcn ipe ipi ipm ipmd ipmnd ipmfg ipmin iput
lags 1 to nlagsrf
end(system)
*
* Determine the maximum estimation range
*
inquire(equation=%modeleqn(reducedform,1)) rstart rend
*
compute ny=%modelsize(reducedform)
*
dec vect[series] factors(nf)
*
* Set up and initialize the D's (one vector per variable)
*
dec vect[vect] d(ny)
do i=1,ny
   compute d(i)=%zeros(nlagsd,1)
end do i
*
* Set up and initialize the C's (one NF x NY matrix per lag)
*
dec vect[rect] c(nlagsc+1)
do i=1,nlagsc+1
   compute c(i)=%zeros(nf,ny)
end do i
*
* Set up and initialize the A's (one NY x NF matrix per lag). This is
* initialized with small random numbers so we'll have derivatives with
* respect to C parameters at the guess values.
*
dec vect[rect] a(nlagsa+1)
do i=1,nlagsa+1
   compute a(i)=%ranmat(ny,nf)*.1
end do i
*
* Put together a full list of variables and lags. For general
* convenience, these are blocked by lag. Also, initialize the (log)
* variable shock variances using a linear regression.
*
dec vect lsv(ny)
dec table fulltable(2,ny*(nlagsrf+1))
do i=1,ny
   compute xseries=%modeldepvars(reducedform)(i)
   linreg(noprint) xseries rstart rend
   # xseries{1 to nlagsd}
   compute lsv(i)=log(%sigmasq)
   compute d(i)=%beta
   do j=0,nlagsrf
      compute fulltable(1,j*ny+i)=xseries
      compute fulltable(2,j*ny+i)=j
   end do j
end do i
equation fulleqn *
# %rlfromtable(fulltable)
*
* These are the multiplied out coefficients in ARF(L)Y(t)=v(t), though
* this is actually organized as ARF(L)Y(t)=-v(t).
*
dec rect arf(ny,ny*(nlagsrf+1))
************************************************
function %%FillARF
local integer i j k f v
compute arf=%zeros(ny,ny*(nlagsrf+1))
do i=1,ny
   compute arf(i,i)=-1.0
   do j=1,nlagsd
      compute arf(i,j*ny+i)=d(i)(j)
   end do j
   do j=0,nlagsa
   	do k=0,nlagsc
         do f=1,nf
            do v=1,ny
               compute arf(i,(j+k)*ny+v)+=a(j+1)(i,f)*c(k+1)(f,v)
            end do v
         end do f
      end do k
   end do j
end do i
end %%FillARF
************************************************
dec symm sv
dec real jacobian
*
* This is called before the main function evaluation to multiply out the
* coefficients and create the target covariance matrix from the log
* variance parameters.
*
function %%IndexStart
compute %%fillarf()
compute sv=%diag(%exp(lsv))
compute jacobian=log(abs(%det(%xsubmat(arf,1,ny,1,ny))))
end %%IndexStart
*
* This includes the normalization that C(L) has a 1 coefficient on lag
* zero for the first variable, and 0 coefficients on the other lags on
* the first variable.
*
nonlin c a d lsv c(1)(1,1)=1.0 c(2)(1,1)=0.0 c(3)(1,1)=0.0
*
* Because the model is linear in variables, non-linear in parameters,
* it's possible to dramatically speed up the calculation by forming the
* crossproduct of the data and computing the implied sigma by a
* quadratic form.
*
cmom(equation=fulleqn) rstart rend
compute rnobs=%nobs
compute %cmom=%cmom/rnobs
*
* Because the likelihood function is computed in one shot using the
* quadratic form, we need to use FIND rather than MAXIMIZE. Since this
* includes current Y's on the RHS, we need to include a Jacobian term,
* which needs to be multiplied by the number of observations, since it
* occurs in each one.
*
find(method=bfgs,iters=800,stderrs) max %logl
   compute %%IndexStart()
   compute [symm] sigmahat=arf*%cmom*tr(arf)
   compute %logl=%logdensitycv(sv,sigmahat,rnobs)+rnobs*jacobian
end find
*
* Create reduced form VAR
*
compute ainv=-1.0*inv(%xsubmat(arf,1,ny,1,ny))
compute rfcoeffs=%zeros(ny*nlagsrf,ny)
do k=1,nlagsrf
   compute phi=ainv*%xsubmat(arf,1,ny,k*ny+1,(k+1)*ny)
   do i=1,ny
      do j=1,ny
         compute rfcoeffs((j-1)*nlagsrf+k,i)=phi(i,j)
      end do j
   end do i
end do k
*
compute %modelsetcoeffs(reducedform,rfcoeffs)
compute factor=ainv*%diag(.5*%exp(lsv))
*
* Compute the decomposition of variance
*
errors(factor=factor,model=reducedform,steps=20)