17/01/2023, 19:14
https://moodle.lse.ac.uk/pluginfile.php/2189838/mod_resource/content/1/HW1.R
library(quantmod)
library(tidyverse)
Copyright By PowCoder代写 加微信 powcoder
library(PerformanceAnalytics)
library(timeSeries)
library(tseries)
library(roll)
library(car)
library(MASS)
library(extraDistr)
library(rugarch)
library(rmgarch)
library(BEKKs)
library(QRM)
library(dplyr)
###########################################################################
# Q1 ———————————————————————-
###########################################################################
# Get stock prices
ENV.CW4 <- new.env() # Create environment where data are stored
Stocks = c('tsla') # Stock names
tickers <- c('tsla') # Stock tickers
# Get the data from YahooFinance
Symbols <- getSymbols(Symbols = tickers, src = 'yahoo',
from = "2013-01-01",
to = "2022-09-30",
env = ENV.CW4)
# Create one XTS object containing adjusted prices of all stocks
Adjusted_Stock_Prices <- do.call(merge, eapply(env = ENV.CW4, Ad))
# Convert into returns ----------------------------------------------------
log_returns <- diff(log(Adjusted_Stock_Prices)) # Comppute daily log returns
log_returns <- na.omit(log_returns) # Remove rows containing na's
AvgRet <- mean(log_returns)
log_returns_demean <- log_returns - AvgRet # Demean returns
###########################################################################
# Q2 ----------------------------------------------------------------------
###########################################################################
https://moodle.lse.ac.uk/pluginfile.php/2189838/mod_resource/content/1/HW1.R 1/4
17/01/2023, 19:14 https://moodle.lse.ac.uk/pluginfile.php/2189838/mod_resource/content/1/HW1.R
par(mfrow=c(1,1))
qqPlot(as.vector(log_returns_demean), xlab = 'Standard Normal Quantiles',
ylab = 'Quantiles of Input Sample', main = names(log_returns_demean),
envelope = FALSE)
# Histogram fitted with normal density
seq_curve <- seq(min(log_returns_demean), max(log_returns_demean), length = 100)
normal_density <- dnorm(x = seq_curve, mean = mean(log_returns_demean), sd = sd(log_returns_demean))
hist(x = log_returns_demean, prob = TRUE, breaks = 100, main = 'Tesla', col = 'blue', xlab = '')
lines(seq_curve, normal_density, lwd = 2, col = 'red')
jarque.bera.test(log_returns_demean)
CV <- qchisq(p = 0.95, df = 2, lower.tail=TRUE)
###########################################################################
# Q3 ----------------------------------------------------------------------
###########################################################################
# Auto correlation function
acf(x = log_returns_demean^2, main = 'Tesla ACF')
lags = c(1, 10, 30)
for (i in 1 : length(lags)) {
print(Box.test(x = log_returns_demean^2, lag = lags[i], type = "Ljung-Box"), fitdf = lags[i])
print(paste('Critical Value is:', qchisq(p = 0.95, df = lags[i], lower.tail=TRUE)))
###########################################################################
# Q4 ----------------------------------------------------------------------
###########################################################################
# EWMA --------------------------------------------------------------------
# Comparing EWMA models with lambda = 0.94
lambda <- c(0.94)
sigmaEWMA <- xts(matrix(nrow = dim(log_returns_demean)[1], ncol = length(lambda)),
order.by = index(log_returns_demean)) # Pre-allocation of estimated conditional volatility matrix
hEWMA <- xts(matrix(nrow = dim(log_returns_demean)[1], ncol = length(lambda)),
order.by = index(log_returns_demean)) # Pre-allocation of estimated conditional variance matrix
hEWMA[1, ] <- var(x = log_returns_demean[1:30]) # Initialize conditional variance by taking sample variance of the first 30 days
par(mfrow=c(1,1))
https://moodle.lse.ac.uk/pluginfile.php/2189838/mod_resource/content/1/HW1.R 2/4
17/01/2023, 19:14 https://moodle.lse.ac.uk/pluginfile.php/2189838/mod_resource/content/1/HW1.R
for (t in 2:dim(log_returns_demean)[1]) {
hEWMA[t] <- lambda * hEWMA[t-1] + (1-lambda) * log_returns_demean[t-1]^2
sigmaEWMA[t] <- sqrt(hEWMA[t])
plot(x = index(sigmaEWMA), y = log_returns_demean, type = 'l', main = paste('EWMA volatility estimate Tesla - lambda =', lambda),
xlab = 'Trading Days', ylab = 'Returns')
lines(x = index(sigmaEWMA), y = 2 * sigmaEWMA, col = 'red')
lines(x = index(sigmaEWMA), y = -2 * sigmaEWMA, col = 'red')
# GARCH(1,1) Model --------------------------------------------------------
# Estimate Garch(1,1) parameters and fit the conditional variance to returns
GARCH_1_1 <- ugarchspec(variance.model = list(model = 'sGARCH', garchOrder = c(1,1)),
mean.model = list(armaOrder = c(0, 0), include.mean = FALSE))
GARCH_1_1_fit <- ugarchfit(spec = GARCH_1_1, data = log_returns_demean, solver = 'hybrid')
hGARCH_1_1 <- xts(x = order.by = index(log_returns_demean))
sigmaGARCH_1_1 <- xts(x = order.by = index(log_returns_demean))
plot(GARCH_1_1_fit) # Running this line shows an interactive menu. Press 1 to plot Series with 2 Conditional SD Superimposed
# GJR-GARCH(1,1) Model --------------------------------------------------------
GJR_GARCH_1_1 <- ugarchspec(variance.model = list(model = 'gjrGARCH', garchOrder = c(1,1)),
mean.model = list(armaOrder = c(0, 0), include.mean = FALSE))
GJR_GARCH_1_1_fit <- ugarchfit(spec = GJR_GARCH_1_1, data = log_returns_demean, solver = 'hybrid')
hGJR_GARCH_1_1 <- xts(x = order.by = index(log_returns_demean))
sigmaGJR_GARCH_1_1 <- xts(x = order.by = index(log_returns_demean))
plot(GJR_GARCH_1_1_fit) # Running this line shows an interactive menu. Press 1 to plot Series with 2 Conditional SD Superimposed
# Log-likelihood test
LR = 2 * (likelihood(GJR_GARCH_1_1_fit) - likelihood(GARCH_1_1_fit))
CV <- qchisq(p = 0.95, df = 1, lower.tail=TRUE)
if (LR > CV) {
print(‘H0: LR = 0 is rejected.’)
print(‘We cannot reject H0’)
###########################################################################
# Q6 ———————————————————————-
###########################################################################
https://moodle.lse.ac.uk/pluginfile.php/2189838/mod_resource/content/1/HW1.R 3/4
17/01/2023, 19:14 https://moodle.lse.ac.uk/pluginfile.php/2189838/mod_resource/content/1/HW1.R
sigmas <- merge(sigmaEWMA, sigmaGARCH_1_1, sigmaGJR_GARCH_1_1) # Merge all volatilities into one object
sigmas_2016 <- sigmas['2016'] # Extract conditional volatilities for 2016
plot(x = index(sigmas_2016), y = log_returns_demean['2016'], type = 'l', main = '3 volatility forecasts for Tesla in 2016', xlab =
'Trading Days', ylab = 'Returns')
lines(x = index(sigmas_2016), y = 2 * sigmas_2016[, 1], col = 'red', lwd = 2)
lines(x = index(sigmas_2016), y = -2 * sigmas_2016[, 1], col = 'red', lwd = 2)
lines(x = index(sigmas_2016), y = 2 * sigmas_2016[, 2], col = 'blue', lwd = 3)
lines(x = index(sigmas_2016), y = -2 * sigmas_2016[, 2], col = 'blue', lwd = 3)
lines(x = index(sigmas_2016), y = 2 * sigmas_2016[, 3], col = 'green', lwd = 1.3)
lines(x = index(sigmas_2016), y = -2 * sigmas_2016[, 3], col = 'green', lwd = 1.3)
abline(v = as.Date('2016-06-22'), lwd = 5, col=alpha(rgb(0.5,0,0.5), 0.3)) # vertical line on 2016-06-22
text(as.Date('2016-07-27'), -.1, as.Date('2016-06-22'))
legend("bottomright", legend = c('EWMA', 'GARCH', 'GJR_GRACH'), lty = 1, col = c("red", "blue", "green"))
https://moodle.lse.ac.uk/pluginfile.php/2189838/mod_resource/content/1/HW1.R 4/4
程序代写 CS代考 加微信: powcoder QQ: 1823890830 Email: powcoder@163.com