new;
library optmum; optset;

@!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!@

@ reading data file @ 

load data[] = c:\gauss50\eco7425\cgrerf.txt; data = reshape(data,rows(data)/3,3);

cg = data[.,1];
re = data[.,2];
rf = data[.,3];

rf = rf/100 + 1;
re = re/100 + 1;

@ computing gross consumption growth rate @

gammac = exp(cg);


n = rows(data);

output file = c:\gauss50\eco7425\cgrerf.out reset ; output on;

@!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!@

declare matrix b;


@!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!@

@ Iterations with different weighting matrices @


count = 1; do while count<=2;

if count == 1;

@ starting parameter values @

b0 = { 0.95, 2};  w = eye(8);


@ Algebraic transformations to handle constraints on parameters @

b0[1] = tan(pi*(b0[1]-0.5)/1);
b0[2] = ln( b0[2] );


elseif count == 2;

b0 = b ;


f = mmcnd(b); 

@ compute the s-hat matrix @

s = f*f' ; s = s/cols(f);

@ computing w matrix as inverse of s-hat @

w = inv(s);


endif;


@!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!@

@ optimisation  @

"optimisation begins"; count;

{b,jmin,grad,retcode} = optmum(&p1,b0);

"optimisation ends"; count;


"count = "; count;
"beta = "; arctan(b[1])*1/pi+0.5;
"sigma = "; exp( b[2] ); 

"jmin = "; jmin;
"gradient = "; grad; 
"retcode = "; retcode;

count = count+1; endo;


@!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!@

@ computing standard errors @

let d[8,2] = 0;   

z=1; do while z<=n-1;
d = d + gradcd(&p2,b);
z=z+1; endo;

d = d/(n-1); 

se = invswp(d'*w*d);  se = sqrt(diag(se));

"standard errors are: "; 
(1/(pi*(1+b[1]^2)))*se[1];
exp( b[2] )*se[2];       

output off;


@!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!@

@ procedure to compute the objective function to minimise @

proc p1(b);
local f, g, jstat;

f = mmcnd(b); 

g = sumc(f')/cols(f);

jstat = g'*w*g;

retp(jstat); endp;


@!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!@

@ procedure to compute derivatives for the covariance matrix @

proc p2(b);
local f;

f = mmcnd(b); 

retp(f[.,z]); endp;


@!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!@

@ procedure to compute the moment conditions @

proc mmcnd(b);
local foce, focf, f;

b[1] = arctan(b[1])*1/pi+0.5;
b[2] = exp( b[2] ) ;


@ Moment Conditions @


@ computing lhs of foc @

foce = b[1]*(gammac^(-b[2])).*re-1;
focf = b[1]*(gammac^(-b[2])).*rf-1;


@ computing the f(i,j) vector @

f = zeros(8,n-1);

f[1,.] = foce[2:n]';
f[2,.] = ( foce[2:n].*gammac[1:n-1] )';
f[3,.] = ( foce[2:n].*re[1:n-1] )';
f[4,.] = ( foce[2:n].*rf[1:n-1] )';
f[5,.] = focf[2:n]';
f[6,.] = ( focf[2:n].*gammac[1:n-1] )';
f[7,.] = ( focf[2:n].*re[1:n-1] )';
f[8,.] = ( focf[2:n].*rf[1:n-1] )';


retp(f); endp;