/* ecn1_2004_2.txt
   program illustrating nonlinear 
   least squares NLS */

closeall;
output file = ecn1_2004_2.out reset;

/* load the data ecn1_wag */

open f1 = ecn1_wag for read;
nr = rowsf(f1);
print "rows of data matrix" nr;
z = readr(f1,nr);

x = ones(nr,1) ~ z[.,2 3 4 5];

y = ln(z[.,1]);

bols = invpd(x'x)*x'y;
print "OLS estimate" bols;


resid = y - x*bols;
r2 = resid.*resid;

proc optobj(delta);
retp( sumc( (r2 - exp(x*delta))^2 ));
endp;


library optmum;
optset;
delta0 = -.1*ones(5,1);

{ deltahat,f,grad,retcode } = optprt(optmum(&optobj,delta0));

/* estimate covariance matrix of disturbances */

varvec = exp(x*deltahat);


/** brute force way to compute fgls 
covmat = diagrv(zeros(nr,nr),varvec);

bfgls = invpd(x'invpd(covmat)*x)*x'invpd(covmat)*y;
**/

/** transfromation method **/
sq_varvec = sqrt(varvec);
ytilda = y ./ sq_varvec ;
xtilda = x ./ sq_varvec;

bfgls = invpd(xtilda'xtilda)*xtilda'ytilda;


print "feasible gls estimates" bfgls;


/* now let us consider the binary dependent variable case discussed
   in class
  one of the problems there was the fact that the linear probability
  model often yielded predicted values outside of the unit interval 

assume that the prob(d=1|X) = exp(Xbeta)/(1+exp(Xbeta))  

as an example, look at the probability of having more than HS
 education given age and Male */

d = z[.,6];
xp = ones(nr,1) ~ z[.,2 3];

fn logit(x) = exp(x) ./ (1+exp(x));

proc logitobj(theta);
local phat,sdev;
phat = logit(xp*theta);
sdev = (d - phat);
retp( sumc( sdev'sdev ) );
endp;

theta0 = 0.1*ones(3,1);

optprt(optmum(&logitobj,theta0));

/* now try with optimally weighted NLS */

proc logitobj_w(theta);
local phat,sdev;
phat = logit(xp*theta);
sdev = (d - phat) ./ sqrt( phat .* (1-phat) );
retp(  sdev'sdev  );
endp;

theta0 = 0.1*ones(3,1);

optprt(optmum(&logitobj_w,theta0));


end;