#************************************************************ 
# GARCH Type Models
# IBM Rate of Returns
# install.packages("tseries")
# install.packages("fGarch")
# install.packages("AER")
# install.packages("fBasics")
#------------------------------------------------------------

library(quantmod)
library(fBasics)
library(TSPred)
library(tseries)
library(fGarch)
library(forecast)
# Yahoo Data in R

getSymbols("IBM",from="2000-01-02", to="2017-06-30")
chartSeries(IBM, theme="white",subset='2007-01::2010-01') 

IBM.rtn = 100*diff(log(IBM$IBM.Adjusted))

chartSeries(IBM.rtn, theme="white",subset='2007::2010')
nob <- nrow(IBM.rtn)
IBM_return <- as.ts(IBM.rtn[2:nob],frequency(365)  )    
head(IBM_return)

# Here is a fund part to estimate various GARCH models 
y <- IBM_return
auto.arima(y, ic = "aic", trace = TRUE)

install.packages("fPortfolio")
library("fPortfolio") 

garchInit=function(series) { 
  Mean=mean(series); Var=var(series); S=1e-6 
  params=c(mu=Mean, omega=0.1*Var, alpha=0.1, beta=0.8) 
  lowerBounds=c(mu=-10*abs(Mean), omega=S^2, alpha=S, beta=S) 
  upperBounds=c(mu=10*abs(Mean), omega=100*Var, alpha=1-S, beta=1-S) 
  cbind(params=params, lowerBounds=lowerBounds, upperBounds=upperBounds) 
} 

init<-garchInit(y) 

# Here is the disrtibution density function for errors
garchDist=function(z,hh) {dnorm(x=z/hh)/hh} 

garchLLH=function(parm, series){ 
  mu=parm[1]; omega=parm[2]; alpha=parm[3]; beta=parm[4] 
  z=(series-mu) 
  Mean=mean(z^2) 
  #Use Filter Representation; 
  e=omega+alpha*c(Mean, z[-length(series)]^2) 
  h=filter(e, beta, "r", init=Mean) 
  hh=sqrt(abs(h)) 
  llh=-sum(log(garchDist(z, hh))) 
  llh 
} 

garchllFit=function(series, params, lowerBounds, upperBounds) { 
  fit=nlminb(start=params, objective=garchLLH, 
             lower=lowerBounds, upper=upperBounds, 
             control=list(trace=3), series=series) 
  fit 
} 

# Here is the nonlinear solutions using PORT routines
garchllFit(y, init[,1], init[,2],init[,3])

# Here is the results from garchFit
g2 <- garchFit(~ garch(1,1), data=y, trace =FALSE)
summary(g2)