/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
/* [ Created with wxMaxima version 0.8.5 ] */

/* [wxMaxima: comment start ]
   Eric Doviak
   original:  22 April 2012
   updated:   30 April 2012

   Math Methods 7025X
   shape of the likelihood function 
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
   First, load the "distrib" package.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
load(distrib)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
   Now, randomly draw 10 values from the standard normal.
   
   y = alpha + beta*x + u
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
N:10$
x:random_normal(0,1,N)$
u:random_normal(0,1,N)$
y:2+3*x+u$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
   Set up the log-likelihood function.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
loglik(alpha,beta,sigma):= -N*log(sigma) - (N/2)*log(2*%pi) - 
(1/2)*sum(((y[i]-alpha-(beta*x[i]))/sigma)^2,i,1,N);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
   Deriviatives must be taken with respect to sigma^2, so define:
   gamma == sigma^2
   and take derivatives with respect to gamma.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
loglik(alpha,beta,gamma):= -(N/2)*log(gamma) - (N/2)*log(2*%pi) - 
(1/(2*gamma))*sum((y[i]-alpha-(beta*x[i]))^2,i,1,N);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
   Maximize it with respect to alpha, beta and gamma. 
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
sol:lbfgs(-loglik(alpha,beta,gamma),'[alpha,beta,gamma],[2.01,2.99,1.01],0.0001,[-1,0])$
alfa_hat:subst(sol[1],alpha)$
beta_hat:subst(sol[2],beta)$
gmma_hat:subst(sol[3],gamma)$

print("")$
print(alpha," = ",alfa_hat)$
print(beta," = ",beta_hat)$
print(sigma^2," = ",gmma_hat)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
   First-order conditions imply that:
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
dla(alpha,beta,gamma):=''(diff(loglik(alpha,beta,gamma),alpha))$
dlb(alpha,beta,gamma):=''(diff(loglik(alpha,beta,gamma),beta))$
dls(alpha,beta,gamma):=''(diff(loglik(alpha,beta,gamma),gamma))$

foc:float(solve([
            dla(alpha,beta,gamma)=0,
            dlb(alpha,beta,gamma)=0,
            dls(alpha,beta,gamma)=0],[alpha,beta,gamma]))$

alfa_foc:subst(foc[1][1],alpha)$
beta_foc:subst(foc[1][2],beta)$
gmma_foc:subst(foc[1][3],gamma)$

print("")$
print(alpha," = ",alfa_foc)$
print(beta," = ",beta_foc)$
print(sigma^2," = ",gmma_foc)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
   Check to see if second-order conditions are satisfied.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
/* set up the Hessian matrix */

daa(alpha,beta,gamma):=''(diff(diff(loglik(alpha,beta,gamma),alpha),alpha))$
dab(alpha,beta,gamma):=''(diff(diff(loglik(alpha,beta,gamma),alpha),beta))$
das(alpha,beta,gamma):=''(diff(diff(loglik(alpha,beta,gamma),alpha),gamma))$

dbb(alpha,beta,gamma):=''(diff(diff(loglik(alpha,beta,gamma),beta),beta))$
dbs(alpha,beta,gamma):=''(diff(diff(loglik(alpha,beta,gamma),beta),gamma))$

dss(alpha,beta,gamma):=''(diff(diff(loglik(alpha,beta,gamma),gamma),gamma))$

H:matrix(
    [daa(alfa_foc,beta_foc,gmma_foc),dab(alfa_foc,beta_foc,gmma_foc),das(alfa_foc,beta_foc,gmma_foc)],
    [dab(alfa_foc,beta_foc,gmma_foc),dbb(alfa_foc,beta_foc,gmma_foc),dbs(alfa_foc,beta_foc,gmma_foc)],
    [das(alfa_foc,beta_foc,gmma_foc),dbs(alfa_foc,beta_foc,gmma_foc),dss(alfa_foc,beta_foc,gmma_foc)])$

print("")$
print("own-partials must be negative:")$
print("")$
print("d^2 loglik"/"(d alpha)^2"," = ",daa(alfa_foc,beta_foc,gmma_foc))$
print("")$
print("d^2 loglik"/"(d beta)^2"," = ",dbb(alfa_foc,beta_foc,gmma_foc))$
print("")$
print("d^2 loglik"/"(d gamma)^2"," = ",dss(alfa_foc,beta_foc,gmma_foc))$
print("")$
print("")$
print("the Hessian matrix:")$
print("H = ",H)$
print("")$
print("determinant of Hessian must be negative")$
print("det(H) = ",determinant(H))$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
   Calculate the standard error of alpha, beta and sigma^2.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
info:-1*invert(H)$

print("")$
print("the information matrix:")$
print("-1*(H^-1) = ",info)$
print("")$
print(alpha," = ",alfa_foc,", se: ",sqrt(info[1,1]))$
print(beta," = ",beta_foc,", se: ",sqrt(info[2,2]))$
print(sigma^2," = ",gmma_foc,", se: ",sqrt(info[3,3]))$
/* [wxMaxima: input   end   ] */

/* Maxima can't load/batch files which end with a comment! */
"Created with wxMaxima"$