*
* @ARMADLM(A=amat,F=fmat,C=Cmat,Z=Zvector) equation
*
* ARMADLM sets up the A and F matrices for a DLM based upon the ARMA
* model in an input equation. This can actually be an ARIMA equation --
* the state-space representation will then have some unit roots.
*
* Parameters:
*
*   equation   The ARMA equation, which you can define with EQUATION
*              or BOXJENK(DEFINE=...)
*
* Options:
*   A=A matrix for input to DLM (output)
*   F=F matrix for input to DLM (output)
*   C=C matrix for input to DLM (output)
*   Z=Z vector for input to DLM (output)
*
* The state space representation is from Jones(1980), "Maximum
* Likelihood Fitting of ARMA Models," Technometrics, pp 389-395. The
* state space representation of x(t) uses a state vector consisting of
* x(t),x(t+1|t), x(t+2|t),...,x(t+r-1|t), where x(s|t) is the best
* linear predictor of x(s) given information through t. For this model,
* A takes the form
*
*  0  1  0 .... 0
*  0  0  1 .... 0
*   ............1
*  p(r-1) p(r-2) ... p(1)
*
* where the p's are the AR coefficients. F is the vector of  "psi
* weights", that is, the responses of the model to a unit shock. The
* dimension r=max(p,q+1). The "C" matrix for this DLM is just
* (1,0,...,0). If there's a CONSTANT in the equation, the "Z" vector is
* (0,...,alpha).
*
* Revision Schedule:
*
*   01/2003 Written by Tom Doan, Estima.
*   04/2008 Renamed as ARMADLM to match the file name.
*   09/2009 Change to produce F matrix, use options for A and F.
*   07/2010 Add C and Z options.
*
procedure ARMADLM eq ap swp
type equation 	eq
type rect 		*ap
type symm 		*swp
*
option rect		*a
option rect		*f
option vect    *c
option vect    *z
*
local rect[int] 	eqtable
local vect 			eqcoeffs
local integer 		p q r i j
local series 		psi
local real        alpha
*
* Scan the equation for the highest AR and MA lags. Also, find the
* coefficient on the CONSTANT, if any.
*
compute alpha=0.0
compute eqtable=%eqntable(eq),eqcoeffs=%eqncoeffs(eq)
compute q=p=0
do i=1,%eqnsize(eq)
   if eqtable(1,i)==%eqndepvar(eq).and.eqtable(2,i)>p
      compute p=eqtable(2,i)
   if eqtable(1,i)==%mvgavge.and.eqtable(2,i)>q
      compute q=eqtable(2,i)
   if eqtable(1,i)==constant
      compute alpha=eqcoeffs(i)
end do i
*
* Figure out the minimal size of the state vector.
*
compute r=%imax(p,q+1)
*
* Compute the impulse responses.
*
impulse(noprint,shock=||1.0||) 1 r
# eq psi
*
* F, C and Z are now known. Set them up if requested by the user.
*
if %defined(f)
   compute f=%xsubvec(psi,1,r)
if %defined(swp)
   compute swp=%outerxx(%xsubvec(psi,1,r))
if %defined(c)
   compute c=%unitv(r,1)
if %defined(z)
   compute z=%unitv(r,r)*alpha
*
* Set up the fixed parts of the A matrix.
*
if %defined(a) {
	dim a(r,r)
	ewise a(i,j)=(j==i+1)
}
if %defined(ap) {
	dim ap(r,r)
	ewise ap(i,j)=(j==i+1)
}
*
* Fill in the AR lags in the bottom row of the A matrix.
*
do i=1,%eqnsize(eq)
   if eqtable(1,i)==%eqndepvar(eq) {
		if %defined(a)
			compute a(r,r+1-eqtable(2,i))=eqcoeffs(i)
		if %defined(ap)
			compute ap(r,r+1-eqtable(2,i))=eqcoeffs(i)
	}
end do i
end

