StatisticalModeling_NASDAQ100.Rmd

---
title: "Financial Time Series - Project"
author: "Shraddha Somani"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Vignette Title}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

## Install Packages
```{r}
library(quantmod)
library(fpp)
library(lattice)
library(timeSeries)
library(rugarch)
library(fUnitRoots)
library(stats)
```

## Get the data
```{r}
getSymbols("^NDX", from="2012-01-01")
NDX[1:10]
dim(NDX)
```

## Chart for NASDAQ100 data
```{r, fig.height=6, fig.width=7.5}
chartSeries(NDX, TA = NULL)
```

## Daily log return
```{r, fig.height=6, fig.width=7}
ndxReturns = dailyReturn(NDX, type = "log")
head(ndxReturns)
length(ndxReturns)

# Remove missing values if any
ndxReturns = na.omit(ndxReturns)
length(ndxReturns)

# ACF and PACF plot for log return
par(mfcol = c(2,1))
acf(ndxReturns)
pacf(ndxReturns)

# Box-test 
Box.test(ndxReturns, type = "Ljung-Box")
```

Analysis:

- From ACF and PACF plot we can see significant trend and seasonality 
- From the box test we can conclude that, the data are not independently distributed; they exhibit serial correlation since the p-value is greater than 0.05

## Create forecast vector
```{r}
winLength = 500
forecastLength = length(ndxReturns) - winLength
forecastLength
forecasts = vector(mode = "character", length = forecastLength)
length(forecasts)
```

## Fit ARIMA model
```{r}
fitarima = auto.arima(ndxReturns)
summary(fitarima)
orderarima = c(length(fitarima$model$phi), 0, length(fitarima$model$theta))
orderarima
```

## Fit GARCH model
```{r}
for (i in 0:forecastLength) {
ndxWinReturns = ndxReturns[(1+i):(winLength + i)]

# Fit GARCH
##--- Method for creating a univariate GARCH specification object prior to fitting. It allows a wide choice in univariate GARCH models , distributions and mean equation modelling---##
##--- Standard specification referred from : https://www.rdocumentation.org/packages/rugarch/versions/1.3-6/topics/ugarchspec-methods ---## 

specification = ugarchspec(variance.model = list(garchOrder=c(1,1)),
                  mean.model = list(armaOrder = c(length(fitarima$model$phi), length(fitarima$model$theta)), include.mean=T),
                  distribution.model="sged")

fitgarch = tryCatch(ugarchfit(specification, ndxWinReturns, solver = 'hybrid'), 
               error=function(e) e, warning=function(w) w)

if(is(fitgarch, "warning")) {
      forecasts[i+1] = paste(index(ndxWinReturns[winLength]), 1, sep="|")
    } else {
      fore = ugarchforecast(fitgarch, n.ahead=1)
      ind = fore@forecast$seriesFor
      forecasts[i+1] = paste(colnames(ind), ifelse(ind[1] < 0, -1, 1), sep="|")
    }    
}

df = data.frame(do.call('rbind', strsplit(as.character(forecasts), '|', fixed = TRUE)))
head(df)
write.csv(df, file="forecasts.csv", row.names = FALSE)        
```

## Convert to time series
```{r}
ndxArimaGarch = xts(df[,-1], order.by = as.Date(df$X1))
class(ndxArimaGarch)
head(ndxArimaGarch)
```

## Create the ARIMA+GARCH returns
```{r, fig.height=6, fig.width=7}
ndxIntersect = merge(ndxArimaGarch[,1], ndxReturns, all = F)
head(ndxIntersect)

ndxArimaGarchReturns = ndxIntersect[,1] * ndxIntersect[,2]
head(ndxArimaGarchReturns)
```

## Backtesting
```{r, fig.height=6, fig.width=7}
ndxArimaGarchCurve = log(cumprod( 1 + ndxArimaGarchReturns))
head(ndxArimaGarchCurve)

ndxBuyHoldCurve = log(cumprod( 1 + ndxIntersect[,2] ) )
head(ndxBuyHoldCurve)

ndxCombinedCurve = merge(ndxArimaGarchCurve, ndxBuyHoldCurve, all = F)
head(ndxCombinedCurve)

xyplot( 
  ndxCombinedCurve,
  superpose=T,
  col=c("darkred", "darkblue"),
  lwd=2,
  key=list( 
    text=list(
      c("ARIMA+GARCH", "Buy & Hold")
    ),
    lines=list(
      lwd=2, col=c("darkred", "darkblue")
    )
  )
)
```