17/01/2023, 19:14
https://moodle.lse.ac.uk/pluginfile.php/2194696/mod_resource/content/1/HW2.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)
# Get data and convert into returns
ENV.CW8 <- new.env() # Create environment where data are stored
Stocks = c('META', 'AMAZON', 'TESLA', 'APPLE', 'NETFLIX', 'ALPHABET') # Stock names
tickers <- c('META', 'AMZN', 'TSLA', 'AAPL', 'NFLX', 'GOOG') # Stock tickers
tickers_cleaned <- c('META', 'AMZN', 'TSLA', 'AAPL', 'NFLX', 'GOOG')
tickers_cleaned <- as.vector(sapply(tickers_cleaned, FUN = function(x) paste(x, '.Adjusted', sep = '')))
# Get the data from YahooFinance
Symbols <- getSymbols(Symbols = tickers, src = 'yahoo',
from = "2013-01-01",
to = "2022-10-15",
env = ENV.CW8)
# Create one XTS object containing adjusted prices of all stocks
Adjusted_Stock_Prices <- do.call(merge, eapply(env = ENV.CW8, Ad))
Adjusted_Stock_Prices <- Adjusted_Stock_Prices[, tickers_cleaned] # Restore the right order of columns
names(Adjusted_Stock_Prices) <- Stocks # Change names from tickers to real names
# Convert into returns
log_returns <- diff(log(Adjusted_Stock_Prices)) # Compute daily log returns
#CalculateReturns(Adjusted_Stock_Prices, method = 'log')
log_returns <- na.omit(log_returns) # Remove rows containing na's
AvgRet <- colMeans(log_returns)
log_returns_demean <- sweep(x = log_returns, MARGIN = 2, STATS = AvgRet) # subtracting column means from each value in each column;
log_returns - AvgRet does not work in R
log_returns_demean
###########################################################################
# Q1 ----------------------------------------------------------------------
https://moodle.lse.ac.uk/pluginfile.php/2194696/mod_resource/content/1/HW2.R 1/3
17/01/2023, 19:14 https://moodle.lse.ac.uk/pluginfile.php/2194696/mod_resource/content/1/HW2.R
###########################################################################
PCA <- prcomp(x = log_returns_demean)
Table_PCA <- rbind(PCA$rotation, summary(PCA)$importance)
par(mfrow=c(1,1))
plot(Table_PCA['Proportion of Variance',], type = 'h', lwd = 5, col = 'blue',
main = 'PC proportions of total variance', xlab = 'PC', ylab = 'Proportion of variance')
###########################################################################
# Q2 ----------------------------------------------------------------------
###########################################################################
PC_loadings <- PCA$x[, 1:nf]
PC_weights <- PCA$rotation[, 1:nf]
spec <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)), # garchOrder = c(q,p) where q = alpha (ARCH) and
p = beta (GARCH)
mean.model = list(armaOrder = c(0, 0), include.mean = FALSE)) # default spec for GARCH(1,1)
mspec <- multispec(replicate(spec, n = nf) )
fitlist <- multifit(multispec = mspec, data = PC_loadings)
sigma_PCs <- sigma(fitlist) # extract info from multifit object using: sigma(), coef(), likelihood(), fitted(), residuals(), show()
var_PCs <- sigma_PCs^2
errors <- log_returns_demean - PC_loadings %*% t(PC_weights)
omega <- diag(x = colMeans(x = errors)^2)
ht <- array(dim = c(length(Stocks), length(Stocks), length(index(log_returns_demean))))
for (i in 1:length(index(log_returns_demean))) {
ht[,,i] <- PC_weights %*% diag(var_PCs[i,]) %*% t(PC_weights) + omega
###########################################################################
# Q3 ----------------------------------------------------------------------
###########################################################################
corr_AAPL_TSLA <- matrix(nrow = length(index(log_returns_demean)), ncol = 1)
for (i in 1:length(index(log_returns_demean))) {
corr_AAPL_TSLA[i] <- ht[3,4,i] / sqrt(ht[3,3,i] * ht[4,4,i])
par(mfrow=c(1,1))
plot(x = index(log_returns_demean), y = corr_AAPL_TSLA, type = 'l', xlab = 'Correlation', ylab = 'Date', main = 'Correlation between
Apple and Tesla')
https://moodle.lse.ac.uk/pluginfile.php/2194696/mod_resource/content/1/HW2.R 2/3
17/01/2023, 19:14 https://moodle.lse.ac.uk/pluginfile.php/2194696/mod_resource/content/1/HW2.R
###########################################################################
# Q4 ----------------------------------------------------------------------
###########################################################################
print(paste('We need', sum(Table_PCA['Cumulative Proportion',] < 0.85) + 1, 'principal components to explain at least 85% of
variance'))
PC_loadings <- PCA$x[, 1:nf]
PC_weights <- PCA$rotation[, 1:nf]
spec <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)), # garchOrder = c(q,p) where q = alpha (ARCH) and
p = beta (GARCH)
mean.model = list(armaOrder = c(0, 0), include.mean = FALSE)) # default spec is GARCH(1,1)
mspec <- multispec(replicate(spec, n = nf) )
fitlist <- multifit(multispec = mspec, data = PC_loadings)
sigma_PCs <- sigma(fitlist) # extract info from multifit object using: sigma(), coef(), likelihood(), fitted(), residuals(), show()
var_PCs <- sigma_PCs^2
errors <- log_returns_demean - PC_loadings %*% t(PC_weights)
omega <- diag(x = colMeans(x = errors)^2)
ht2 <- array(dim = c(length(Stocks), length(Stocks), length(index(log_returns_demean))))
for (i in 1:length(index(log_returns_demean))) {
ht2[,,i] <- PC_weights %*% diag(var_PCs[i,]) %*% t(PC_weights) + omega
###########################################################################
# Q5 ----------------------------------------------------------------------
###########################################################################
cov_AAPL_TSLA <- matrix(nrow = length(index(log_returns_demean)), ncol = 2)
for (i in 1:length(index(log_returns_demean))) {
cov_AAPL_TSLA[i, 1] <- ht[3,4,i]
cov_AAPL_TSLA[i, 2] <- ht2[3,4,i]
plot(x = index(log_returns_demean), y = cov_AAPL_TSLA[,1], type = 'l', col = 'blue', xlab = 'Date', ylab = 'Covariance', main =
'Covariances of Apple and Tesla', lwd = 4)
lines(x = index(log_returns_demean), y = cov_AAPL_TSLA[,2], col = 'red')
legend("topleft", legend = c('Cov - 6 PCs', 'Cov - 3 PCs'), lty = 1, col = c("blue", "red"))
https://moodle.lse.ac.uk/pluginfile.php/2194696/mod_resource/content/1/HW2.R 3/3
程序代写 CS代考 加微信: powcoder QQ: 1823890830 Email: powcoder@163.com