*
* @DLMGLS( options ) depvar start end
* #  list of regressors
*
* DLMGLS estimates a linear equation by GLS where the error process follows a
* dynamic linear model with structure similar to those which are supported by the
* RATS instruction DLM.
*
* The (univariate) error process e(t) takes the form
*
*  X(t) = A X(t-1) + W(t)
*  e(t) = c'X(t) + v(t)
*
* Options (all must be constants or fixed matrices):
*   A = A matrix from state space representation of error process
*   C = c vector from state space representation of error process
*   SW = SYMMETRIC covariance matrix of W
*   SX0 = SYMMETRIC pre-sample covariance matrix of X
*   SV  = variance of v [0.0]
*
* Revision Schedule:
*   05/2007 Written by Tom Doan, Estima
*   10/2008 Several corrections made. Added PRINT/NOPRINT option.
*
procedure DLMGLS depvar start end
type series depvar
type integer start end
*
option rect   A
option vect   C
option symm   SW
option symm   sx0
option real   sv        0.0
option switch print     1
option switch vweighted 1
*
local vect[vect] x
local integer    n m startl endl
local integer    time j
local equation   eqn
local rect[int]  regtable
local symm       cmom sx
local real       csc
local vect       sc
local vect       z
local vect       xy beta
local symm       xx
local real       rss seesq
*
compute n=%rows(A)
*
* Get the explanatory variables for the regression into eqn
*
equation eqn depvar
compute m=%eqnsize(eqn)
compute regtable=%eqntable(eqn)
*
inquire(equation=eqn) startl<<start endl<<end
*
dim x(m+1) z(m+1)
*
compute cmom=%zeros(m+1,m+1)
do j=1,m+1
   compute x(j)=%zeros(n,1)
end do j
compute sx=sx0
do time=startl,endl
   if %defined(A)
      compute sx=A*sx*tr(A)
   if %defined(SW)
      compute sx=sx+sw
   compute csc=%qform(sx,c)+sv
   compute sc =sx*c
   compute sx =sx-(1.0/csc)*%outerxx(sc)
   do j=1,m+1
      if %defined(A)
         compute x(j)=A*x(j)
      if j==m+1
         compute z(j)=depvar(time)-%dot(c,x(j))
      else
      if regtable(1,j)==0
         compute z(j)=1.0-%dot(c,x(j))
      else
         compute z(j)=([series] regtable(1,j))(time-regtable(2,j))-%dot(c,x(j))
      compute x(j)=x(j)+(z(j)/csc)*sc
      if vweighted
         compute z(j)=z(j)/sqrt(csc)
   end do j
   disp "At time" time z
   compute cmom=cmom+z*tr(z)
end do time
disp "Components" cmom
compute xx   =inv(%xsubmat(cmom,1,m,1,m))
compute xy   =%xsubmat(cmom,1,m,m+1,m+1)
compute beta =xx*xy
compute rss  =cmom(m+1,m+1)-%dot(%beta,xy)
compute seesq=rss/(%nobs-m)
compute xx   =xx*seesq
*
linreg(create,print=print,equation=eqn,covmat=xx,noscale,coeffs=beta,$
   title="Regression with State Space Errors") * startl endl
end







