objects();                     # prints a list of all variables in memory
remove(list=objects());        # sets objects() as a list, and then removes this list
options(scipen=20);           # sets the decimal places in scientific notation
sink(file = "arch.out", append = FALSE, type = "output", split = TRUE);
# creates a disk file of the terminal/screen output; type can be changed to
# log inputted commands as well.


ARCH.1 <- function(theta,x)
{
logl = 0;
e    <<- 0;                          # variable made global
sigma2 <<- 0;
npar <<- 3;
n    <<- nrow(x);
mu     = theta[1];
omega  = theta[2];
alpha  = theta[3];
uncond = omega/(1-alpha);            # unconditional volality


i = 1;
while (i <= n){
  e[i] <<- x[i,1] -mu;
  if (i <= 1){
    sigma2[i] <<- uncond;
  }else{                                  
    sigma2[i] <<- omega + alpha*(e[i-1]^2);
  }
  if(sigma2[i] <= 0.0){
    return(10.0^200);                   #penalize non-positive variances
  }
  
  logl = logl -.5*log(2*pi) -.5*log(sigma2[i]) -.5*(e[i]^2)/sigma2[i];
  i = i + 1;
}

return(-logl);
}


# MAIN PROGRAM

myData <- read.table("archdata.txt", header=TRUE);
y      = as.matrix(myData$perc100);



print("NOBS=");
print(nrow(y));

result <- nlm(ARCH.1, c(.001,.001,.001), y, hessian=T, steptol=.0000001, stepmax=100,
              iterlim=100,print.level=2,
              fscale=500,typsize=c(1.0,.1,.1));

# note there are a few optimizers available for R.  I prefer nlm() over optim().
# Also note that the optimizer defaults to the BFGS algoritm and minimizes the
# negative of the log-likelihood.

print(result);


invhess <- solve( result$hessian );
print( cbind(result$estimate,sqrt( diag(invhess) )) );


par(mfrow=c(2,1))         # multifigure setup: 2 rows, 1 cols
ts.plot(y,main="Returns");
ts.plot(sigma2,main="Conditional Variance");
par(ask = TRUE);

## Some density plot for standardized residuals

z <- e/sqrt(sigma2);

par(mfrow=c(2,1));        # set up the graphics 
d <- density(z,n=2000);          #kernel plot of z
plot(d,main="Kernel density of z and implied normal density");
normx <- seq(min(z),max(z),length=50);
normy <- dnorm(normx,mean=mean(z),sd=sd(z));
lines(normx,normy,col="blue");
qqnorm(z);             # normal Q-Q plot  
qqline(z);             # add a line    
par(ask = TRUE);


d$y <- log(d$y);
normy <- log(normy);
plot(d,main="LOG Kernel density of z and implied normal density");
#normx <- seq(min(z),max(z),length=50);
#normy <- dnorm(normx,mean=mean(z),sd=sd(z));
lines(normx,normy,col="blue");

## perform LB tests
cat("\nLB tests on z and z^2\n");
print( Box.test(z, lag = 1, type = c("Ljung-Box")) );
print( Box.test(z, lag = 10, type = c("Ljung-Box")) );
z2 <- z^2;
print( Box.test(z2, lag = 1, type = c("Ljung-Box")) );
print( Box.test(z2, lag = 10, type = c("Ljung-Box")) );

cat("\nLB tests on returns^2\n");
print( Box.test(y^2, lag = 1, type = c("Ljung-Box")) );


plot(lag(e),sqrt(sigma2));
par(ask = TRUE);


# students may wish to compare results to the garch() R canned routine included
# in the t-series package.

#sink(file = NULL)