*
* @VARDLM(options) LagMats
*
* Generates the state-space system matrices for a vector autoregressive
* process whose lag coefficients are in the parameter LagMats
* VECTOR[RECT]. You can also use it to set up the framework (proper
* dimensions, all the fixed values put in place) based upon an input
* number of lags (LAGS option) to fill out the non-fixed elements later.
*
* The simplest state-space (dynamic linear model) representation of
*
*  Y(t)=A1 Y(t-1) + A2 Y(t-2) + ... + Ap Y(t-p) + U(t)
*
* where the Y's and U's are N-vectors and the A's are N x N lag matrices
*
* has state-vector X(t)=[Y(t),Y(t-1),...,Y(t-p+1] (thus p x N total states, blocked by lags)
*
*     |A1,A2,...,Ap-1,Ap|
*     |I , 0,...,0   ,0 |
* A = |0 , I,...,0   ,0 |
*     | ...     ,0   ,0 |
*     |0 , 0,...,I   ,0 |
*
* where all the elements are N x N matrix blocks
*
* F = (transpose of) |I,0,0,...,0| Np x N matrix consisting of p blocks of N x N matrices.
*
* C = same as F
*
* Parameters:
*
*  LagMats (optional, input)  VECT[RECT] of VAR coefficient matrices on consecutive lags
*
*  Options:
*
*  LAGS=number of lags in model (input)
*    Can be used as an alternative to the LagMats parameter to set up the
*    framework of the matrices if you don't yet know the coefficients.
*
*  N=number of variables in model (input)
*    Can be used if LagMats is empty (for a zero lag VAR) to give the size
*    of the set of variables.
*
*  A=DLM system "A" matrix (output; input/output if used with POKE)
*  F=DLM system "F" matrix (output)
*  C=DLM system "C" matrix (output)
*
*  POKE/[NOPOKE]
*    If POKE, assumes that the A option matrix has the framework set up and
*    only changes the first row to insert the coefficients out of LagMats.
*    This avoids repeating the set up of the fixed elements.
*
* Revision Schedule:
*  07/2025 Written by Tom Doan, Estima
*
procedure VARDLM LagMats
type vect[rect] LagMats
option rect *A
option rect *F
option rect *C
option int  lags 0
option int  n    0
option switch poke 0
*
local int inputLags workLags l i j nw adim states
*
if %defined(LagMats)
   compute inputLags=%size(LagMats)
else
   compute inputLags=0
*
compute workLags=%imax(inputLags,lags)
*
* If lags are zero, we need A to be an N x N zero matrix.
*
if inputLags>0
   compute nw=%rows(LagMats(1))
else
   compute nw=n
*
* This is the number of lags in the state vector. This has to be at
* least one.
*
compute stateLags=%imax(workLags,1)
compute states=nw*stateLags
*
if %defined(A) {
   if .not.poke
      compute A=%zeros(states,states)
   do l=1,workLags
      do i=1,nw
         do j=1,nw
            if l<=%size(LagMats)
               compute A(i,(l-1)*nw+j)=LagMats(l)(i,j)
            if l<>1.and.i==j
               compute A((l-1)*nw+i,(l-2)*nw+j)=1
         end do j
      end do i
   end do l
}
if poke
   return

if %defined(F) {
   dim F(states,nw)
   do l=1,stateLags
      do i=1,nw
         do j=1,nw
            compute F((l-1)*nw+i,j)=(l==1.and.i==j)
         end do j
      end do i
   end do l
}

if %defined(C) {
   dim C(states,nw)
   do l=1,stateLags
      do i=1,nw
         do j=1,nw
            compute C((l-1)*nw+i,j)=(l==1.and.i==j)
         end do j
      end do i
   end do l
}
end
