Research project - R programming / analytics

Assignment 3

Sukhman Kaur

2021-06-09

"Part A – Risks 1. Using two securities from your portfolio (from different sectors)"

## [1] "Part A – Risks\n1. Using two securities from your portfolio (from different sectors)"

"a) Construct a variance-covariance matrix:"

## [1] "a) Construct a variance-covariance matrix:"

r = getOption("repos") r["CRAN"] = "http://cran.us.r-project.org" options(repos = r)

library(quantmod)

## Loading required package: xts

## Loading required package: zoo

## ## Attaching package: 'zoo'

## The following objects are masked from 'package:base': ## ## as.Date, as.Date.numeric

## Loading required package: TTR

## Registered S3 method overwritten by 'quantmod': ## method from ## as.zoo.data.frame zoo

library(xts) getwd() # verify your working directory

## [1] "/Users/sukhmankaur/Downloads"

# Importing Amazon data from Yahoo Finance data.AMZN.2 = read.csv("~/Downloads/AMZN-2.csv", stringsAsFactors=TRUE)

# Check data head(data.AMZN.2)

## Date Open High Low Close Adj.Close Volume ## 1 2010-06-09 120.31 121.47 117.36 117.91 117.91 7369200 ## 2 2010-06-10 120.00 123.50 119.20 123.21 123.21 6061800 ## 3 2010-06-11 121.39 123.53 120.29 123.03 123.03 4204600 ## 4 2010-06-14 124.24 125.70 123.50 123.83 123.83 3923000 ## 5 2010-06-15 123.20 126.92 122.50 126.84 126.84 4541000 ## 6 2010-06-16 125.39 127.98 125.36 126.90 126.90 3964300

# Variable Date date = as.Date(data.AMZN.2$Date,format="%M-%D-%Y")

# Replace date variable with the date data.AMZN.2 = cbind(date,data.AMZN.2[,-1])

# Sort data in chronological order data.AMZN.2 = data.AMZN.2[order(data.AMZN.2$date),] class(data.AMZN.2)

## [1] "data.frame"

# Renaming the columns names(data.AMZN.2) = paste(c("AMZN.2.Open", "AMZN.2.High", "AMZN.2.Low", "AMZN.2.Close", "AMZN.2.Adjusted", "AMZN.2.Volume")) data.AMZN.2[c(2:5,nrow(data.AMZN.2)),]

## AMZN.2.Open AMZN.2.High AMZN.2.Low AMZN.2.Close AMZN.2.Adjusted ## 2 <NA> 120.00 123.50 119.20 123.21 ## 3 <NA> 121.39 123.53 120.29 123.03 ## 4 <NA> 124.24 125.70 123.50 123.83 ## 5 <NA> 123.20 126.92 122.50 126.84 ## 2517 <NA> 2500.20 2530.00 2487.34 2524.06 ## AMZN.2.Volume NA ## 2 123.21 6061800 ## 3 123.03 4204600 ## 4 123.83 3923000 ## 5 126.84 4541000 ## 2517 2524.06 3970700

"Since we want forth-year data, in my case it’s the year 2009."

## [1] "Since we want forth-year data, in my case it’s the year 2009."

BTC.USD <- read.csv("~/Downloads/BTC-USD.csv")

# create variable date date = as.Date(BTC.USD$Date,format="%Y-%m-%d")

# Replace date variable with the date data.BTC.USD = cbind(date,BTC.USD[,-1])

# Sort data in chronological order data.BTC.USD = data.BTC.USD[order(data.BTC.USD$date),] class(data.BTC.USD)

## [1] "data.frame"

# Renaming the columns names(data.BTC.USD) <- paste(c("BTC.USD.Open", "BTC.USD.High", "BTC.USD.Low", "BTC.USD.Close", "BTC.USD.Adjusted", "BTC.USD.Volume")) data.BTC.USD[c(1:5,nrow(data.BTC.USD)),]

## BTC.USD.Open BTC.USD.High BTC.USD.Low BTC.USD.Close BTC.USD.Adjusted ## 1 2015-06-09 228.537994 230.953995 227.929001 229.048004 ## 2 2015-06-10 228.994995 229.781998 228.009995 228.802994 ## 3 2015-06-11 228.854996 230.287003 228.766998 229.705002 ## 4 2015-06-12 229.705002 231.057007 229.313004 229.981995 ## 5 2015-06-13 229.919998 232.651993 229.210007 232.401993 ## 1828 2020-06-09 9774.360352 9836.369141 9664.719727 9795.700195 ## BTC.USD.Volume NA ## 1 229.048004 28353100 ## 2 228.802994 15904800 ## 3 229.705002 14416000 ## 4 229.981995 14017700 ## 5 232.401993 13305300 ## 1828 9795.700195 23717842783

# Using the subset command for getting data. test.BTC.USD <-subset(BTC.USD, + index(BTC.USD) >= "2015-07-01" &+ index(BTC.USD) <= "2016-07-31") test.BTC.USD[c(1:3,nrow(test.BTC.USD))]

## [1] Date Open High ## <0 rows> (or 0-length row.names)

##BTC.USD.ret = Delt(BTC.USD$BTC.USD.Adjusted) ##BTC.USD.ret[c(1:3,nrow(BTC.USD.ret)),] "Combine the two return series"

## [1] "Combine the two return series"

##"Create a vector of weights AMZN: 25% ##WGT.2asset <- c(0.25,0.75) ##WGT.2asset <- matrix(WGT.2asset,1) ##WGT.2asset

“Part B – Suitable distributions for returns”

“1. Fit the data using GHD, HYP and NIG” #plot the density function - BITCOIN

library(fBasics)

## Loading required package: timeDate

## Loading required package: timeSeries

## ## Attaching package: 'timeSeries'

## The following object is masked from 'package:zoo': ## ## time<-

## ## Attaching package: 'fBasics'

## The following object is masked from 'package:TTR': ## ## volatility

library(timeSeries)

##Port.ret<-0.25*AMZN.2.ret +0.75*BTC.USD.ret ##Port.ret<-Port.ret[-1,] ##ef = density(Port.ret, na.rm = TRUE) ##par(mfrow = c(1,1)) ##plot(ef)

#Fit the Generalized Hyperbolic Distribution

library(ghyp)

## Loading required package: numDeriv

## Loading required package: MASS

#Fit the Hyperbolic Distribution

hypfit<- fit.hypuv(Port.ret, symmetric = FALSE, control = list(maxit = 1000), na.rm = TRUE)

#Fit the Normal Inverse Gaussian Distribution

nigfit <- fit.NIGuv(Port.ret, symmetric =FALSE, control = list(maxit = 1000), na.rm= TRUE)

“2. Plot the combined density functions”

ghddens <- dghyp(efx, hypfit) nigdens <- dghyp(efx, mean = mean(Port.ret, na.rm=TRUE), sd = sd(c(Port.ret[, 1]), na.rm=TRUE)) col.def <- c(“black”, “red”, “blue”, “green”, “orange”) plot(ef, xlab = "“, ylab = expression(f(x)), ylim = c(0, 30)) lines(ef$x, ghddens, col = "red") lines(ef$x, hypdens, col =”blue“) lines(ef$x, nigdens, col = "green") lines(ef$x, nordens, col =”orange“) legend(”topleft“, legend = c(”empirical“,”GHD“,”HYP“,”NIG“,”NORM"), col = col.def, lty = 1)

“3. Create a Q-Q plot”

qqghyp(ghdfit, line = TRUE, ghyp.col = “red”, plot.legend = FALSE, gaussian = FALSE, main = "“, cex = 0.8) qqghyp(hypfit, add = TRUE, ghyp.pch = 2, ghyp.col =”blue“, gaussian = FALSE, line = FALSE, cex = 0.8) qqghyp(nigfit, add = TRUE, ghyp.pch = 3, ghyp.col =”green“, gaussian = FALSE, line = FALSE, cex = 0.8) legend(”topleft“, legend = c(”GHD“,”HYP“,”NIG"), col = col.def[-c(1,5)], pch = 1:3)

“4. Make a model recommendation using lik.ratio.test”

AIC <- stepAIC.ghyp(Port.ret,control = list(maxit = 1000))

head(AIC$fit.table)

“Use the function “lik.ratio.test” to perform a likelihood-ratio test on fitted generalized hyperbolic distribution objects of class mle.ghyp."

“The likelihood-ratio test can be used to check whether a special case of the generalized hyperbolic distribution is the “true” underlying distribution."

“statistic: the value of the L-statistic.”

“p.value: the p-value for the test (the p-value is less than 0.05, then there is evidence against the null hypothesis)”

“df: the degrees of freedom for the L-statistic”

“H0: a boolean stating whether the null hypothesis is TRUE or FALSE (if TRUE there is no relationship between the data sets)”

LRghdnig <- lik.ratio.test(ghdfit, nigfit) LRghdnig

LRghdhyp <- lik.ratio.test(ghdfit, hypfit) LRghdhyp

“5. Calculate and plot the VaR (using all models)”

p <- seq(0.001, 0.05, 0.001)

ghd.VaR <- abs(qghyp(p,ghdfit)) hyp.VaR <- abs(qghyp(p, hypfit)) nig.VaR <- abs(qghyp(p, nigfit)) nor.VaR <- abs(qnorm(p, mean = mean(Port.ret, na.rm=TRUE), sd = sd(c(Port.ret[, 1]), na.rm = TRUE))) emp.VaR <- abs(quantile(x = Port.ret,probs= p, na.rm = TRUE)) plot(emp.VaR, type = “l”, xlab = "“, ylab =”VaR“, axes = FALSE, ylim = range(c(hyp.VaR, nig.VaR, ghd.VaR, nor.VaR, emp.VaR))) box() axis(1, at = seq(along = p), labels = names(emp.VaR), tick = FALSE) axis(2, at = pretty(range(emp.VaR, ghd.VaR, hyp.VaR, nig.VaR, nor.VaR))) lines(seq(along = p), ghd.VaR, col =”red“) lines(seq(along = p), hyp.VaR, col =”blue“) lines(seq(along = p), nig.VaR, col =”green“) lines(seq(along = p), nor.VaR, col =”orange“) legend(”topright“, legend = c(”Empirical“,”GHD“,”HYP“,”NIG“,”Normal"), col = col.def, lty = 1)

“6. Calculate and plot the ES (using all models)”

ghd.ES <- abs(ESghyp(p,ghdfit)) hyp.ES <- abs(ESghyp(p, hypfit)) nig.ES <- abs(ESghyp(p, nigfit)) nor.ES <- abs(mean(Port.ret, na.rm=TRUE) - sd(c(Port.ret[, 1]), na.rm = TRUE) dnorm(qnorm(1 - p)) / p) obs.p <- ceiling(plength(Port.ret)) emp.ES <- sapply(obs.p, function(x) abs(mean(sort(c(Port.ret)) [1:x]))) plot(emp.ES, type = “l”, xlab = "“, ylab =”ES“, axes = FALSE, ylim = range(c(hyp.ES, nig.ES, ghd.ES, nor.ES, emp.ES), na.rm = TRUE)) box() axis(1, at = 1:length(p), labels = names(emp.VaR), tick = FALSE) axis(2, at = pretty(range(emp.ES, ghd.ES, hyp.ES, nig.ES, nor.ES))) lines(1:length(p), ghd.ES, col =”red“) lines(1:length(p), hyp.ES, col =”blue“) lines(1:length(p), nig.ES, col =”green“) lines(1:length(p), nor.ES, col =”orange“) legend(”topright“, legend = c(”Empirical“,”GHD“,”HYP“,”NIG“,”Normal"), col = col.def, lty = 1,cex = 0.7)