#==============
# Preliminaries
#==============

    # -------------------------------------------
    # Nuclear power station example - Example 5.2
    # from Brazzale, Davison and Reid (2007)
    #
    # This file does the OLS fit
    # -------------------------------------------
library(marg)
data(nuclear)

require(MASS)
require(nlme)
require(numDeriv)
source("hessian.txt")



#===============================
 loglik <- function( beta, y, x) 
#===============================
{
        # Log-likelihood function for nuclear power station 
        # data - OLS model

    mu = beta[1] + beta[2]*x[,1] + beta[3]*x[,2] + 
         beta[4]*x[,3] + beta[5]*x[,4] + beta[6]*x[,5] + 
         beta[7]*x[,6] 
    sum( (y-mu)*(y-mu) / 2)
}


    # ----------------------
    # Least squares fit in R
    # ----------------------
nuc.norm <- lm( log(cost) ~ date + log(cap) + NE + CT + log(N) + PT,
                data=nuclear)
summary(nuc.norm)

attach(nuclear)
y = log(cost)
x = cbind(date,log(cap),NE,CT,log(N),PT)
detach(nuclear)

param = nuc.norm$coef

truese = hess2corr(ginv(vcov(nuc.norm) / summary(nuc.norm)$sigma^2))

cc <- chkerrors(loglik, param, y, x, truese)
print(cc)