*
* @DSGEControl( options ) input *output
*
* Procedure for doing optimal control on a linear rational expectations
* model. Algorithm from Dennis(2007), "Optimal Policy in Rational
* Expectations Models: New Solution Algorithms", Macroeconomic Dynamics,
* vol 11, 31-55.

* Parameters:
*   input = input components generated by DSGE with form=second. This puts
*   the model into the form
*
*    I(1)X(t)=I(2)X(t-1)+I(3)EX(t+1)+I(4)w(t)+I(5)
*
*   output = output components after optimal control in the same form as
*   the inputs.
*
* Options:
*   BETA=discount rate in the cost functional
*   NX  =number of controls[1]. These should be the final variables in the
*     variables list on DSGE.
*   Q   =quadratic form matrix on the controls [zero matrix]
*   R   =quadratic form matrix on the controlled variables [zero matrix]
*
procedure DSGEControl comps revised
type vect[rect] comps *revised
*
option integer nx  	1
option real    beta  .99
*
option rect    q
option rect    r
*
local integer ny ix iy il in
local rect    a0 a1 a2
local rect    b0 b1 b2
local rect    a5 a6
*
compute ny=%rows(comps(1))-nx
compute a0=%xsubmat(comps(1),1,ny,1,ny)
compute a1=%xsubmat(comps(2),1,ny,1,ny)
compute a2=%xsubmat(comps(3),1,ny,1,ny)
compute b0=-1.0*%xsubmat(comps(1),1,ny,ny+1,ny+nx)
compute b1=%xsubmat(comps(2),1,ny,ny+1,ny+nx)
compute b2=%xsubmat(comps(3),1,ny,ny+1,ny+nx)
compute a5=%xsubmat(comps(4),1,ny,1,%cols(comps(4)))
compute a6=%xsubmat(comps(5),1,ny,1,1)
*
compute ix=1,iy=ix+nx,il=iy+ny,in=il-1+ny
*
dim revised(5)
compute revised(1)=%zeros(in,in)
compute %psubmat(revised(1),il,iy,a0)
compute %psubmat(revised(1),il,ix,-1.0*b0)
compute %psubmat(revised(1),iy,il,tr(a0))
if %defined(r)
   compute %psubmat(revised(1),iy,iy,r)
compute %psubmat(revised(1),ix,il,-1.0*tr(b0))
if %defined(q)
   compute %psubmat(revised(1),ix,ix,q)
*
compute revised(2)=%zeros(in,in)
compute %psubmat(revised(2),il,iy,a1)
compute %psubmat(revised(2),iy,il,(1.0/beta)*tr(a2))
compute %psubmat(revised(2),ix,il,(1.0/beta)*tr(b2))
compute %psubmat(revised(2),il,ix,b1)
*
compute revised(3)=%zeros(in,in)
compute %psubmat(revised(3),ix,il,beta*tr(b1))
compute %psubmat(revised(3),il,iy,a2)
compute %psubmat(revised(3),il,ix,b2)
compute %psubmat(revised(3),iy,il,beta*tr(a1))
*
compute revised(4)=%zeros(in,%cols(a5))
compute %psubmat(revised(4),il,1,a5)
*
compute revised(5)=%zeros(in,1)
compute %psubvec(revised(5),il,a6)
end


