#load Rmetrics packages
library(mvtnorm)
library(mgcv)
library(nnet)
library(MASS)
library(fBasics)
library(fCalendar)
library(fSeries)
library(fEcofin)
library(fMultivar)
library(fPortfolio)
library(fOptions)
library(fExtremes)

#imports data from CSV file
stock.data = read.csv("stock_data.csv", header=T)

#extracts vector of dates and vector of returns
dates = stock.data[,3]
brk.returns = stock.data[,5]
spx.returns = stock.data[,6]

#create time series and plot
brkSeries = timeSeries(data=brk.returns, charvec=dates)
spxSeries = timeSeries(data=spx.returns, charvec=dates)

plot(brkSeries)
plot(spxSeries)

#plot Ljung Box, ACF, and Standardized Residuals
tsdiag(arima(brkSeries, c(0,0,0)), gof.lag=20)
tsdiag(arima(spxSeries, c(0,0,0)), gof.lag=20)

#estimate ARMA model to return series
brk.arima = arima(brkSeries, c(0, 0, 1))
spx.arima = arima(spxSeries, c(0, 0, 1))

tsdiag(spx.arima, gof.lag=20)

#form squared return series and plot
brk.var = brk.returns*brk.returns
spx.var = spx.returns*spx.returns

brkVarSeries = timeSeries(data=brk.var, charvec=dates)
spxVarSeries = timeSeries(data=spx.var, charvec=dates)

plot(brkVarSeries)
plot(spxVarSeries)

#plot Ljung Box, ACF, and Standardized Residuals for squared returns
tsdiag(arima(brkVarSeries, c(0,0,0)), gof.lag=20)
pacf(brkVarSeries)

tsdiag(arima(spxVarSeries, c(0,0,0)), gof.lag=20)
pacf(spxVarSeries)

#fit ARMA model to variance series
brkvar.arima = arima(brkVarSeries, c(1, 0, 0))
tsdiag(brkvar.arima, gof.lag=20)

spxvar.arima = arima(spxVarSeries, c(7, 0, 0))
tsdiag(spxvar.arima, gof.lag=20)

#test for normality
qqgaussPlot(brkSeries)
qqgaussPlot(spxSeries)
jarqueberaTest(brkSeries)
jarqueberaTest(spxSeries)

#calculate skewness and kurtosis
skewness(brk.returns)
skewness(spx.returns)
kurtosis(brk.returns)
kurtosis(spx.returns)