*
*  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.
*
/* Below is a sample program demonstrating how to use the procedure:

cal 1947 1 12
all 1988:3
open data rates.rat
data(format=rats)
*
dec rect a(1,1)
dec vect c(1)
dec symm sx0(1,1) sw(1,1)
compute rho=.75
*
compute a=||.75||
compute c=||1.0||
compute sx0=||1.0/(1-.75**2)||
compute sw=||1.0||
*
source(noecho) dlmgls.src
@dlmgls(a=a,c=c,sx0=sx0,sw=sw) fygl 1948:1 1988:3
# constant fygm3{0 to 3}
*/

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
*
local vect[vect] x
local integer n m
local integer time j
local equation eqn
local rect[int] regtable
local symm cmom sx
local real csc
local vect sc
local vect z
*
compute n=%rows(A)
equation eqn depvar
compute m=%eqnsize(eqn)
compute regtable=%eqntable(eqn)
*
dim x(m+1) z(m+1) cmom(m+1,m+1)
*
compute cmom=%const(0.0)
do j=1,m+1
   dim x(j)(n)
   compute x(j)=%const(0.0)
end do j
compute sx=sx0
do time=start,end
   if %defined(A)
      compute sx=A*sx*tr(A)
   if %defined(SW)
      compute sx=sx+sw
   disp sx
   compute csc=%qform(sx,c)+sv
   compute sc =sx*tr(c)
   compute sx =sx-(1.0/csc)*%outerxx(sc)
   disp csc sc sx
   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
      compute z(j)=z(j)/sqrt(csc)
   end do j
   compute cmom=cmom+z*tr(z)
end do time
compute beta=inv(%xsubmat(cmom,1,m,1,m))*%xsubmat(cmom,1,m,m+1,m+1)
end

