R codes

R.codes

###############################
## Access the relevant files ##
###############################

returns <- read.csv("Returns_CNX_500.csv")
returns1 <- returns
nifty <- read.csv("Nifty_returns.csv")
mibor <- read.csv("MIBOR.csv", na.strings="#N/A")
exchange <- read.csv("Exchange_rates.csv", na.strings="#N/A")

###################################################################
## Dealing with missing values ##
###################################################################

for(i in 2:ncol(returns))
{
  returns1[, i] <- approx(returns$Year, returns1[ ,i], returns$Year)$y  # Replacing the missing values by linear approximates
 
}

## Dealing with blanks in the MIBOR rates ##

mibor[, 2] <- approx(as.Date(mibor$Dates, '%d-%b-%y'), mibor[ ,2], as.Date(mibor$Dates, '%d-%b-%y'))$y
for(k in 2:nrow(mibor))
{
  mibor$Change1[k] <- ((mibor$MIBOR[k] - mibor$MIBOR[k-1])/mibor$MIBOR[k-1])*100
}

## Dealing with blanks in the exchange rates ##

exchange[, 2] <- approx(as.Date(exchange$Year,'%d-%b-%y'), exchange[ ,2], as.Date(exchange$Year, '%d-%b-%y'))$y
exchange$Change <- as.numeric(exchange$Change)
for(j in 2:nrow(exchange))
{
exchange$Change[j] <- ((exchange$Exchange.rates[j] - exchange$Exchange.rates[j-1])/exchange$Exchange.rates[j-1])*100
}


#######################################
### Summary Statistics calculations ###
#######################################

library(moments)
library(tseries)


## Unit root testing ##


adf.test(exchange$Exchange.rates)
adf.test(exchange$Change)
adf.test(mibor$MIBOR)
adf.test(mibor$Change1)
adf.test(nifty$S...P.Cnx.Nifty)

skewness(mibor$MIBOR)
skewness(mibor$Change1)
skewness(exchange$Exchange.rates)
skewness(exchange$Change)
skewness(nifty$S...P.Cnx.Nifty)


kurtosis(mibor$MIBOR)
kurtosis(mibor$Change1)
kurtosis(exchange$Exchange.rates)
kurtosis(exchange$Change)
kurtosis(nifty$S...P.Cnx.Nifty)

var(exchange$Exchange.rates)
var(exchange$Change)
var(mibor$MIBOR)
var(mibor$Change1)
var(nifty$S...P.Cnx.Nifty)

Box.test(exchange$Exchange.rates, type = "Box-Pierce", lag= 10)
Box.test(exchange$Change, type = "Ljung-Box", lag= 10)
Box.test(mibor$MIBOR, type = "Ljung-Box", lag= 10)
Box.test(mibor$Change1, type = "Ljung-Box", lag= 10)
Box.test(nifty$S...P.Cnx.Nifty, type = "Ljung-Box", lag= 10)


#####################################################
## Calculating the principal component of returns ##
#####################################################

# Store the relevent data as a matrix and then pass it through princomp()

ret <- as.matrix(returns1[ ,-1], nrow = dim(returns1)[1], ncol = dim(returns1)[2]-1)
princ.return <- princomp(ret)
plot(princ.return, main = "Principal component variance") # Plot of variance of the components
barplot(height=princ.return$sdev[1:10]/princ.return$sdev[1], main = "Relative variance of the principal components" ) # Relative variance of the components 
total.var <- sum(princ.return$sdev^2)
pc1.var <- princ.return$sdev[1]^2

pc10.var <- sum(princ.return$sdev[1:10]^2)

## Loadings of first 10 components ##

# Extracting the components using the loadings of the PCA

load.cp <- loadings(princ.return)[,1:10] 
pr.cp <- ret %*% load.cp
pr <- as.data.frame(pr.cp)

##########################################
## Regressions single factor model CAPM ##
##########################################


reg.capm <- lapply( 1:10, function (i) { lm(pr[, i] ~ nifty$S...P.Cnx.Nifty) })

# Storing the coefficients of the 10 regressions in "cfs"

cfs <- sapply(reg.capm, function(x) as.vector(x$coefficients))
format(cfs, digits=3)

# Similarly storing the R^2 of the 10 regressions in 

r.sq.values <- lapply(reg.capm, function(x) summary(x)$r.squared)
format(r.sq.values, digits = 3)


colnames(cfs) <- lapply(1:10, function(x) paste("Comp", as.character(x), sep=''))
rownames(cfs) <- c("Intercept", "Nifty")

# Storing the output in a LaTeX friendly format

write.table(format(cfs, digits=3), row.names=T, col.names=T, sep=" & ", eol=" \\\\ \n", file="my.temp.file.txt")
write.table(format(r.sq.values, digits = 3), file="my.temp.file.txt", append=T, col.names=F, row.names=T, sep = " & ", eol="\\\\ \n")


###################################
## Regressions multi factor APT ##
###################################


reg.apt <- lapply( 1:10, function (i) { lm(pr[, i] ~ nifty$S...P.Cnx.Nifty + mibor$Change1 + exchange$Change) })
tmp.apt <- reg.apt[[1]]

# Storing the coefficients and the R^2 values

cfs1 <- sapply(reg.apt, function(x) as.vector(x$coefficients))
r.sq.values1 <- lapply(reg.apt, function(x) summary(x)$r.squared)

colnames(cfs1) <- lapply(1:10, function(x) paste("Comp", as.character(x), sep=''))
rownames(cfs1) <- c("Intercept", "Nifty", "MIBOR" , "INR/USD Exchange rate")

write.table(format(cfs1, digits=4), row.names=T, col.names=T, sep=" & ", eol=" \\\\ \n", file="my.temp.file.APT.txt")
write.table(format(r.sq.values1, digits = 4), file="my.temp.file.APT.txt", append=T, col.names=F, row.names=T, sep = " & ", eol="\\\\ \n")

#############################
## Plotting various series ##
#############################

pr$dates <- as.data.frame(mibor$Dates)

## Non stationary independent variables ##

pdf("indep_var_ns.pdf", height = 7.5, width = 10)
par(mfrow = c(2, 1))
plot(as.Date(mibor$Dates,'%d-%b-%y'), mibor$MIBOR, xlab= "Date", 
     ylab= "3-month MIBOR rates (%age)", type='l', col='red', 
     main="3-month MIBOR rates")
abline(h = 0, lty = 8, col = "gray")
plot(as.Date(exchange$Year, '%d-%b-%y'), exchange$Exchange.rates, xlab= "Date", 
     ylab= "IND/USD Exchange rates", type='l', col='red', 
     main="IND/USD Exchange rate")
abline(h = 0, lty = 8, col = "gray")
dev.off()

## Stationary independent variables ##

pdf("indep_var.pdf", height = 7.5, width = 10)
par(mfrow = c(3, 1))
plot(as.Date(mibor$Dates,'%d-%b-%y'), mibor$Change1, xlab= "Date", 
     ylab= "Change in 3-month MIBOR(%age)", type='l', col='royalblue', 
     main="%age change in MIBOR rates")
abline(h = 0, lty = 8, col = "gray")
plot(as.Date(nifty$Date,'%d-%b-%y'), nifty$S...P.Cnx.Nifty, xlab= "Date", 
     ylab= "NIFTY returns(%age)", type='l', col='royalblue', 
     main="NIFTY returns")
abline(h = 0, lty = 8, col = "gray")
plot(as.Date(exchange$Year, '%d-%b-%y'), exchange$Change, xlab= "Date", 
     ylab= "IND/USD Exchange rates change(%age)", type='l', col='royalblue', 
     main="IND/USD Exchange rate changes(%age)")
abline(h = 0, lty = 8, col = "gray")
dev.off()


plot(as.Date(exchange$Year[seq(from=1 , to=length(exchange$Year), by=5)],'%d-%b-%y'), exchange$Change[seq(from=1 , to=length(exchange$Change), by=5)], xlab= "Date", ylab= "Exchange rates change", type='l', col='red', lwd=1)
tmp <- seq(from=1, to=length(..), by=5)
plot(as.Date(e$Y[tmp]), e$E[tmp])


# Plot of the dependents
max <- max(pr$Comp.1)/11
plot(as.Date(mibor$Dates,'%d-%b-%y')[seq(from=1 , to=length(mibor$Dates), by=10)], (-pr$Comp.1/max)[seq(from=1 , to=length(pr$Comp.1), by=10)] , xlab= "Date", 
     ylab= "PC 1 / NIFTY", type='l', col='red', 
     main="First Principal component and NIFTY")
lines(as.Date(mibor$Dates,'%d-%b-%y')[seq(from=1 , to=length(mibor$Dates), by=10)], nifty$S...P.Cnx.Nifty[seq(from=1 , to=length(mibor$Dates), by=10)], col = 'royalblue', type = "l")

cor(pr$Comp.1, nifty$S...P.Cnx.Nifty   )

plot(as.Date(mibor$Dates,'%d-%b-%y')[seq(from=1 , to=length(mibor$Dates), by=10)], nifty$S...P.Cnx.Nifty[seq(from=1 , to=length(mibor$Dates), by=10)], col = 'royalblue', type = "l",
     ylab = "Nifty returns", xlab = "Date", main = "Nifty returns" )

###################################################
#### Comparision of the regression F-Statistic ####
###################################################


# ANOVA table for illustration from where we can extract the F-values for the second
# regression directly

tmp1 <- reg.capm[[1]]; tmp2 <- reg.apt[[1]]
anova(tmp1, tmp2)

all.fstats <- mapply(anova, reg.capm, reg.apt )
lapply(all.fstats, function(x) x$F )
all.fstats[[6]]

# Extracting the F-stats from the ANOVA table of all the 10 regressions

all.Fs <- mapply(function(x, y) anova(x,y)$F[2], reg.capm, reg.apt)
all.Ps <- mapply(function(x, y) anova(x,y)[[ 2, "Pr(>F)"]], reg.capm, reg.apt)

names(anova(tmp1, tmp2))


##################################################################################
######################## TIME SERIES TERM PAPER ##################################
##################################################################################

###### Plotting the residuals after regression ####

install.packages("lmtest")
library(lmtest)

# Extract the residuals from all the regressions

res.apt <- sapply(reg.apt, function(x) as.vector(x$residuals))
res.capm <- sapply(reg.capm, function(x) as.vector(x$residuals))


colnames(res.apt) <- lapply(1:10 , function(x) paste("Residuals reg", as.character(x) , sep= ' '))
colnames(res.capm) <- lapply(1:10, function(x) paste("Residuals reg", as.character(x) , sep= ' '))

# Perform the Lyung-box test on the residuals of all the regressions both CAPM and APT

lm.res <- function(res){ Box.test(res, lag = 10, type = "Ljung-Box")}
box.res <- lapply(1:10, function(x) lm.res(res.apt[, x]))

lm.res.capm <- function(res){ Box.test(res, lag = 10, type = "Ljung-Box")}
box.res.capm <- lapply(1:10, function(x) lm.res.capm(res.capm[, x]))

# Extracting the values of the Q-Stat and the corrosponding p-value

stats <- lapply(box.res, function(x) x$statistic)
p.vals <- lapply(box.res, function(x) x$p.value)

stats.capm <- lapply(box.res.capm, function(x) x$statistic)
p.vals.capm <- lapply(box.res.capm, function(x) x$p.value)


#################################################################################
###################### MLE FOR GARCH(1,1)########################################
#################################################################################

# Getting all the data in one data frame

data.garch <- data.frame(Component1 = pr$Comp.1, Exchange.change = exchange$Change, Niftychange = nifty$S...P.Cnx.Nifty, mibor.change = mibor$Change1)
write.csv(data, file = "GARCH-Data", col.names = TRUE)



# Specifying functions:

CalcResiduals <- function(th, data) {
   # Calculates the e_t and h_t for the GARCH(1, 1) model with given parameters.
   #
   # Args:
   #   th: Parameters
   #          th[1] -> mean
   #          th[2] -> alpha.0
   #          th[3] -> alpha.1
   #          th[4] -> beta.1
   #          th[5] -> sigma.0
   #          th[6] -> a
   #          th[7] -> b
   #          th[8] -> c
   #   data: The input data
   #
   # Returns: A list containing et and ht.
   
 th[1] -> mean
 th[2] -> alpha.0
 th[3] -> alpha.1
 th[4] -> beta.1
 th[5] -> sigma.0
 th[6] -> a
 th[7] -> b
 th[8] -> c

 
 if(is.null(data$Niftychange) || is.null(data$Exchange.change) || is.null(data$mibor.change)) print("nifty column not present")

# These are the residuals "y"

y = data$Component1 - a*data$Niftychange - b*data$Exchange.change - c*data$mibor.change


 n <- length(y)
 sigma.sqs <- vector(length=n)
 sigma.sqs[1] <- sigma.0 ^ 2
 for(ii in c(1:(n-1))) { ## This loop is where the h_t are calculated
   sigma.sqs[ii + 1] <- (
     alpha.0 +
     alpha.1 * (y[ii] - mean) ^ 2 + 
     beta.1 * sigma.sqs[ii])
 }

 return(list(et = y, ht = sigma.sqs)) # Returns the list of e_t and h_t
}



GarchLogL <- function(th, data) {
 # Calculates the Log-Likelihood of the GARCH(1, 1) model with given
 # parameters. It is intended for use with nlp or other optimization
 # routines to arrive at the best GARCH model. This can also be called for
 # subsets of the data and then summed to arrive at other models.
 #
 # Args:
 #   th : This is a vector containing the parameters for the model:
 #     th[1] -> mean
 #     th[2] -> alpha.0
 #     th[3] -> alpha.1
 #     th[4] -> beta.1
 #     th[5] -> sigma.0
 #     th[6] -> a
 #     th[7] -> b
 #     th[8] -> c
 # Returns:
 #   The negative conditional log likelihood of the model
 

 res <- CalcResiduals(th, data)
 sigma.sqs <- res$ht
 y <- res$et

# Assuming normal density of the errors dnorm() gives the density of normal dist.
# this will return the negative of the log likelihood.
 
 return (-sum(dnorm(y[-1], mean=th[1] , sd=sqrt(sigma.sqs[-1]), log=TRUE))) 
}



## Main program where we will call all our previous written functions

GarchLogLSimpl <- function(th, y) { # Only setting the mean of the errors to 0  
    GarchLogL(c(0, th), y)
  }

GarchLogLSimpl2 <- function(th, y) { # Setting mean and b,c = 0
    GarchLogL(c(0, th, 0,0), y)
  }

GarchLogLSimpl3 <- function(th, y) { # Setting mean of errors, alpha.1, beta.1 = 0
    GarchLogL(c(0, th[1], 0,  0,1, th[2], th[3], th[4]), y)
  }

# The nlm() function performs the non-linear optimization with "p" as the initial 
# values of the parameters and then goes ahead with the iterations till the value
# of the likelihood function does not converge.

system.time(fit2 <- nlm(GarchLogLSimpl, p = rep(1,7), hessian = TRUE, data <- data.garch , iterlim = 500, print.level = 2 ))
sqrt(diag(solve(fit2$hessian)))

# Squre root of the Diagonal of the inverse of the hessian matrix gives the standard
# errors of the estimates


system.time(fit3 <- nlm(GarchLogLSimpl2, p = rep(1,5), hessian = TRUE, data <- data.garch , iterlim = 500, print.level = 2 ))
sqrt(diag(solve(fit2$hessian)))

system.time(fit4 <- nlm(GarchLogLSimpl3, p = rep(1,4), hessian = TRUE, y = data.garch , iterlim = 500, print.level = 2 ))
sqrt(diag(solve(fit4$hessian)))

# Checking for remaining ARCH effect in the standardized shocks

aaa <- CalcResiduals(c(0, 5.9560857,   0.1524731,   0.8143006,  13.5448251, -12.2535545,   1.4845980,
 0.2730385) , data = data.garch)

std.shock.1 <- aaa$et/(aaa$ht)^0.5 
Box.test(std.shock.1, lag = 10, type = "Ljung-Box")

good