/*
** Lesson 10.7: Maximum Likelihood Estimation
** AR(1), MA(1), ARMA(1,1)
*/
use gpe2;
output file = gpe\output10.7 reset;
load data[40,6]= gpe\cjx.txt;
     
year = data[2:40,1];
X = ln(data[2:40,2]);
L = ln(data[2:40,3]);
K = ln(data[2:40,5]);
data=X~L~K;

@ OLS estimates as initial values @
b=data[.,1]/(ones(rows(data),1)~data[.,2:3]);

call reset;
_nlopt=1;
_method=0;
_iter=100;
_conv=1;

_b=b|0.5;
_jacob=&jcb;
_names = {"CONSTANT","LN(L)","LN(K)","AR(1)"};
call estimate(&ar,data);

_b=b|0;
_jacob=0;
_names = {"CONSTANT","LN(L)","LN(K)","MA(1)"};
call estimate(&ma,data);

_b=b|0.5|0;
_jacob=&jcb;
_names = {"CONSTANT","LN(L)","LN(K)","AR(1)","MA(1)"};
call estimate(&arma,data);

end;

proc jcb(x,b);  @ jacobian for AR(1) and ARMA(1,1) @
    local j;
    j=ones(rows(x),1);
    j[1]=sqrt(1-b[4]^2);
    retp(j);
endp;

proc ar(x,b);
    local n,e,u;
    n=rows(x);
    e=x[.,1]-b[1]-b[2]*x[.,2]-b[3]*x[.,3];
    u=e-b[4]*lagn(e,1);
    @ first obs transformation @
    u[1]=sqrt(1-b[4]^2)*e[1];
    retp(u);
endp;

proc ma(x,b);
    local n,e,u;
    n=rows(x);
    e=x[.,1]-b[1]-b[2]*x[.,2]-b[3]*x[.,3];
    u=recserar(e,e[1],b[4]); @ 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[4]*u[i-1];
        i=i+1;
    endo;
*/
    retp(u);
endp;

proc arma(x,b);
    local n,e,u,v;
    n=rows(x);
    e=x[.,1]-b[1]-b[2]*x[.,2]-b[3]*x[.,3];
    u=e-b[4]*lagn(e,1);
    @ first obs transformation @
    u[1]=sqrt(1-b[4]^2)*e[1];
    v=recserar(u,u[1],b[5]);
    retp(v);
endp;