/*
** ARMA1.GPE - nonlinear maximum likelihood estimation
**             for AR(1), MA(1), and ARMA(1,1) models
**
** (C) Copyright 2001 by Applied Data Associates
** All Rights Reserved.
**
** This Software Product is PROPRIETARY SOURCE CODE OF APPLIED
** DATA ASSOCIATES. This File Header must accompany all files
** using any portion, in whole or in part, of this Source Code.
** This Software Product is designed to be used with GPE2 and
** GAUSS. If you use this Source Code for research and development,
** a proper reference is required. If you wish to distribute any 
** portion of the proprietary Source Code, in whole or in part, you 
** must first obtain written permission from Applied Data Associates.
**
** Consider the linear regression equation y = xb + e, and e is
** assumed to be one of three basic autoregressive structures: AR(1),
** MA(1), ARMA(1,1). The predefined residual functions ar1, ma1, and 
** arma1 are given for these three models respectively. To estimate a 
** linear regression model with AR(1) error structure, just call 
** estimate procedure with ar1 residual function as follows:
** 
** ==> call estimate(&ar1,z);
**
** Similarly, to estimate MA(1) model,
**
** ==> call estimate(&ma1,z);
**
** To estimate ARMA(1,1) model,
**
** ==> call estimate(&arma1,z);
**
** The data matrix z =[y,x] is arranged so that the left-hand-side 
** dependent variable of the linear regression equation y = xb + e
** is in the first column, followed by the right-hand-side explanatory
** variables including contant. Furthermore, to estimate AR(1) and
** ARMA(1,1) models, the Jacobian function must be defined for the
** control variable _jacob as:
**
** ==> _jacob=&jf1;
**
** This is done before calling estimate procedure with the relevant
** ar1 or arma1 function. For MA(1) model, the Jacobian is not needed.
**
** The predefined first-order autoregressive models are for the linear 
** regression equation only. It is straightforward to extend the model 
** for nonlinear regression equation and for higher order of ARMA 
** strucutre. The application module ARMA1.GPE is designed that it can
** be easily modified and extended.
**
*/ 

proc jf1(x,b);   @ jacobian for AR(1) and ARMA(1,1) @
    local j,n,k;
    n=rows(x);   @ number of observations @
    k=cols(x)-1; @ number of regressors @
    j=ones(n,1); @ make sure AR(1) parameter is b[k+1]  @
    j[1]=sqrt(1-b[k+1]^2);
    retp(j);
endp;

proc ar1(x,b);
    local n,k,e,u;
    n=rows(x);
    k=cols(x)-1;                @ rows(b)=k+1 @
    e=x[.,1]-x[.,2:k+1]*b[1:k]; @ b[k+1] is AR(1) parameter @
    u=e-b[k+1]*lagn(e,1);
    u[1]=sqrt(1-b[k+1]^2)*e[1]; @ first obs transformation @
    retp(u);
endp;

proc ma1(x,b);
    local n,k,e,u;
    n=rows(x);
    k=cols(x)-1;                @ rows(b)=k+1 @
    e=x[.,1]-x[.,2:k+1]*b[1:k]; @ b[k+1] is MA(1) parameter @
    u=recserar(e,e[1],b[k+1]);  @ u[1]=e[1] since u[0]=0 @
/*
    @ recursive computation of errors using @
    @ built-in RECSERAR is the same as below: @
    u=e;    @ initialize: u[1]=e[1] @
    i=2;
    do until i>n;
        u[i]=e[i]+b[k+1]*u[i-1];
        i=i+1;
    endo;
*/
    retp(u);
endp;

proc arma1(x,b);
    local n,k,e,u,v;
    n=rows(x);
    k=cols(x)-1;                @ rows(b)=k+2 @
    e=x[.,1]-x[.,2:k+1]*b[1:k];
    u=e-b[k+1]*lagn(e,1);
    @ first obs transformation @ 
    u[1]=sqrt(1-b[k+1]^2)*e[1]; @ b[k+1] is AR(1) parameter @
    v=recserar(u,u[1],b[k+2]);  @ b[k+2] is MA(1) parameter @
    retp(v);
endp;