Predictive Analysis of Stock Market Performance Using Statistical and Financial Methods
2022-12-03
Section 0 - Synopsis
The aim of this project is to analyze and compare the daily closing prices of ten selected stocks and the S&P 500 index during two different periods, with the objective of predicting the next day’s return, Value at Risk (VaR), and Expected Shortfall (ES) for each stock. The two periods are January 1, 2018, to November 30, 2022, and January 1, 2013, to November 30, 2017 respectively.
The analysis employs statistical and financial theorems including ordinary least square analysis and principal component regression analysis to build the model and exponential smoothing to increase accuracy in capturing volatility. The effectiveness of the model was evaluated through the calculation of Sharpe ratios. Additionally, risk management techniques were employed to assess potential risks.
Section 1 - Introduction
Background
The Standard and Poor’s 500, or simply the S&P 500, is a stock market index tracking the stock performance of 500 large companies listed on stock exchanges in the United States.
This study tracks and investigates the predicted next-day S&P500 return, which can help potential investors or risk managers to make investments.
Objective and Data Description
The objective of this project is to predict the next day’s S&P 500 return by utilizing a lagged S&P 500 and the returns of the daily closing prices of ten selected stocks over the past five years. The data was obtained from the Yahoo Finance source via the quantmod R package. The ten companies, namely Apple Inc, Microsoft Corp, Amazon.com Inc, Alphabet Inc, Tesla Inc, NVIDIA Corp, PepsiCo Inc, Meta Platforms Inc, Costco Wholesale Corp, and Broadcom Inc, were selected as they constitute the top 100 constituents by weight and collectively account for approximately 20% of the total market capitalization. The data used for analysis spans from January 1, 2018, to November 30, 2022, and encompasses the daily closing prices and S&P 500 index price for the past five years.
In this study, a modified prediction algorithm was developed and the dataset was partitioned into three distinct sets: 50% for training, 25% for validation, and 25% for testing to achieve the research objective.
The read.bossa.data, pred.footsie.prepare,
first.acf.squares.train, vol.exp.sm, three
thresh.reg, Sharpe.curves functions were
obtained from the lecture notes of a financial statistics course (code
ST326) at the London School of Economics and Political Science. They
were authored by Professor Clifford Lam, and have undergone several
modifications. For an in-depth understanding of the functions, please do
not hesitate to contact me, and I will furnish you with a comprehensive
explanation.
Load packages:
library(quantmod)
library(tseries)
library(rugarch)
library(forecast)
library(ggplot2)Section 2 - EXPLORATORY DATA ANALYSIS
Comparison of two economic period
Due to the global impact of the COVID-19 pandemic over the past five years, the world economy has been impacted to varying degrees. Two sets of data were compared throughout this report, contrasting the data of the same companies during the period of 2013 to 2017, when the economy was in an upward trend, to the present situation.
Importing data
getSymbols(c("AAPL","MSFT","AMZN","GOOG","AVGO","TSLA","NVDA","PEP", "META","COST","^GSPC"), src = "yahoo", from = as.Date("2013-01-01"), to = as.Date("2017-12-01"))## [1] "AAPL" "MSFT" "AMZN" "GOOG" "AVGO" "TSLA" "NVDA" "PEP" "META"
## [10] "COST" "^GSPC"AAPL2 <- AAPL
MSFT2 <- MSFT
AMZN2 <- AMZN
GOOG2 <- GOOG
AVGO2 <- AVGO
TSLA2 <- TSLA
NVDA2 <- NVDA
PEP2 <- PEP
META2 <- META
COST2 <- COST
GSPC2 <- GSPC
getSymbols(c("AAPL","MSFT","AMZN","GOOG","AVGO","TSLA","NVDA","PEP", "META","COST","^GSPC"), src = "yahoo", from = as.Date("2018-01-01"), to = as.Date("2022-12-01"))## [1] "AAPL" "MSFT" "AMZN" "GOOG" "AVGO" "TSLA" "NVDA" "PEP" "META"
## [10] "COST" "^GSPC"str(AAPL)## An 'xts' object on 2018-01-02/2022-11-30 containing:
## Data: num [1:1238, 1:6] 42.5 43.1 43.1 43.4 43.6 ...
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : chr [1:6] "AAPL.Open" "AAPL.High" "AAPL.Low" "AAPL.Close" ...
## Indexed by objects of class: [Date] TZ: UTC
## xts Attributes:
## List of 2
## $ src : chr "yahoo"
## $ updated: POSIXct[1:1], format: "2023-03-06 00:37:50"Checking for missing values
#justify for missing values
any(is.na(AAPL[,4]))## [1] FALSEany(is.na(AMZN[,4]))## [1] FALSEany(is.na(AVGO[,4]))## [1] FALSEany(is.na(COST[,4]))## [1] FALSEany(is.na(GOOG[,4]))## [1] FALSEany(is.na(GSPC[,4]))## [1] FALSEany(is.na(META[,4]))## [1] FALSEany(is.na(MSFT[,4]))## [1] FALSEany(is.na(NVDA[,4]))## [1] FALSEany(is.na(PEP[,4]))## [1] FALSEany(is.na(TSLA[,4]))## [1] FALSEany(is.na(AAPL2[,4]))## [1] FALSEany(is.na(AMZN2[,4]))## [1] FALSEany(is.na(AVGO2[,4]))## [1] FALSEany(is.na(COST2[,4]))## [1] FALSEany(is.na(GOOG2[,4]))## [1] FALSEany(is.na(GSPC2[,4]))## [1] FALSEany(is.na(META2[,4]))## [1] FALSEany(is.na(MSFT2[,4]))## [1] FALSEany(is.na(NVDA2[,4]))## [1] FALSEany(is.na(PEP2[,4]))## [1] FALSEany(is.na(TSLA2[,4]))## [1] FALSE#so there are no missing valuesThe is.na() function proves no missing value in the
document.
Export data into text files-1
Export data into text files-2
The study employed the quantmod package to obtain stock market data,
which was then prepared for analysis using the
read.bossa.data function. This function reads stock market
data from files in the specified directory and organizes the data into
matrices. A vector of names of the stock files to be read is taken as an
input to the function.
The read.bossa.data function was used to create a log-price plot of the stocks, which was subsequently exported into text files for further analysis.
Log-price plot of 2018-2022
source("myfunctions.R")
###plot the log-price graph
vec.names = c("AAPL","MSFT","AMZN","GOOG","AVGO","TSLA","NVDA","PEP","META","COST","GSPC")
ind <- read.bossa.data(vec.names)
mn <- max(ind$starting.indices)
mx <- dim(ind$arranged.closes)[2]
ts.plot(t(ind$arranged.closes[1:8,mn:mx]), col=1:8)
par(new=TRUE)
ts.plot(t(ind$arranged.closes[9:11,mn:mx]), lty=3, col=1:3)
legend("bottomleft", vec.names, lty=c(rep(1,8),rep(3,3)), col=1:11)
In this study, the trend of the daily closing prices of the selected ten companies was analyzed. With the exception of META, all companies displayed an upward trend, with Tesla exhibiting exceptional growth. The COVID-19 pandemic in 2019 had a varying impact on the prices of these companies, with Tesla being the most affected and Costco being minimally impacted. In 2022, a slight economic downturn was observed, resulting in declining prices for most companies, except for Costco, which showed an increase.
Log-price plot of 2013-2017
###plot the log-price graph
vec.names = c("AAPL2","MSFT2","AMZN2","GOOG2","AVGO2","TSLA2","NVDA2","PEP2","META2","COST2","GSPC2")
ind <- read.bossa.data(vec.names)
mn <- max(ind$starting.indices)
mx <- dim(ind$arranged.closes)[2]
ts.plot(t(ind$arranged.closes[1:8,mn:mx]), col=1:8)
par(new=TRUE)
ts.plot(t(ind$arranged.closes[9:11,mn:mx]), lty=3, col=1:3)
legend("bottomleft", vec.names, lty=c(rep(1,8),rep(3,3)), col=1:11)
Based on the historical closing price plot of these stocks from
2013-2017, we can observe varying trends and patterns across different
stocks. For instance, some stocks such as Broadcom Inc (AVGO) and Nvidia
(NVDA) experienced significant growth and appreciation over this period.
Others, such as Pepsi (PEP) and Costco (COST), had more stable and
consistent growth. The S&P 500 index (^GSPC) has also shown growth
over this period, while Tesla (TSLA) and Apple (AAPL) experienced
significant volatility and fluctuations in their stock prices.
Overall, these plots can provide valuable insights into the performance and volatility of these stocks over the given period, which can be attributed to a variety of factors such as changes in market demand, company announcements, economic indicators, government policies, and global events. By analyzing the patterns and trends in the plot, investors and analysts can gain a better understanding of the underlying forces driving the performance of these stocks and make informed investment decisions.
Section 3 - Volatility study
pred.footsie.prepare For the 11 daily return series, I write an R programme to use exponential smoothing to estimate their daily volatilities over the horizon of the training data. Individual λs for each series are estimated by MLE.
The function pred.footsie.prepare prepares the data for
the prediction exercise and splits them into a 50% train, 25% validation
and 25% test sets. The function takes three input parameters:
max.lag, split, and mask.
max.lag specifies the maximum lag to be included in the
prediction. split indicates how much of the data to include
in the training and validation sets, and mask specifies
which other indices to include (1 for yes, 0 for no).
shift.indices represents whether the closing values of each
stock are on the same day (1 for previous day, 0 for same day).
The function pred.footsie.prepare2 is further
transferred for data from 2013-2017.
In this report, it should be noted that all selected companies have their closing values on the same day. Therefore, the function is defined with all shift indices set to 0, meaning that data up to t, which is the current time, is being considered.
Function of volatility
Concept for using exponential smoothing to estimate daily volatilities
Set GARCH(1,1) model for each variables
The observation that volatility is correlated has led to the study of GARCH (1,1) models.
For the individual lambdas:
Here I built a GARCH (1,1) model for each time series. I define lp as
the log-price for each company and differencing to eliminate the first
observation with lp. r = r[2:length(r)] is used to remove
the effect of NAs produced in the previous step, and then use the
garch() function in tseries package to define the
distribution.
In the r code, fit$coe[3] is used to output the value of
λ, which is exactly the coefficient of \(\sigma^2_{t-1}\) estimated by MLE.
Here is an example of Apple(AAPL), and the detailed code could be find in the appendix.
#Apple Inc
lp.AAPL = log(AAPL[,5])
r1 = diff(lp.AAPL)
r1 = r1[2:length(r1)]
fit = garch(r1)##
## ***** ESTIMATION WITH ANALYTICAL GRADIENT *****
##
##
## I INITIAL X(I) D(I)
##
## 1 7.539065e-02 1.000e+00
## 2 5.000000e-02 1.000e+00
## 3 5.000000e-02 1.000e+00
##
## IT NF F RELDF PRELDF RELDX STPPAR D*STEP NPRELDF
## 0 1 -9.348e+02
## 1 4 -9.359e+02 1.16e-03 4.78e-03 6.1e-02 4.5e+04 1.0e-02 1.07e+02
## 2 5 -9.387e+02 2.94e-03 3.46e-03 6.1e-02 2.0e+00 1.0e-02 5.16e+01
## 3 6 -9.418e+02 3.38e-03 4.44e-03 1.3e-01 2.0e+00 2.0e-02 3.46e+01
## 4 7 -9.460e+02 4.38e-03 9.14e-03 1.9e-01 2.0e+00 4.0e-02 9.71e+00
## 5 8 -9.486e+02 2.77e-03 7.08e-03 1.3e-01 2.0e+00 4.0e-02 1.15e+00
## 6 9 -9.510e+02 2.48e-03 3.64e-03 1.1e-01 2.0e+00 4.0e-02 1.34e-01
## 7 10 -9.517e+02 7.92e-04 8.60e-04 1.1e-01 2.0e+00 4.0e-02 6.23e-02
## 8 11 -9.518e+02 6.58e-05 1.13e-03 2.0e-01 1.9e+00 8.0e-02 1.83e-02
## 9 12 -9.523e+02 5.05e-04 1.26e-03 8.9e-02 1.9e+00 4.0e-02 4.69e-03
## 10 13 -9.524e+02 1.53e-04 3.58e-04 8.6e-02 1.9e+00 4.0e-02 1.41e-03
## 11 14 -9.525e+02 5.17e-05 1.43e-04 8.7e-02 1.7e+00 4.0e-02 1.13e-03
## 12 15 -9.525e+02 7.06e-06 5.84e-05 2.9e-02 0.0e+00 1.4e-02 5.84e-05
## 13 16 -9.525e+02 2.47e-05 2.99e-05 2.9e-02 8.1e-02 1.4e-02 3.00e-05
## 14 18 -9.525e+02 2.06e-06 3.14e-06 4.4e-03 7.3e-01 2.1e-03 4.29e-06
## 15 19 -9.525e+02 1.35e-07 1.54e-07 2.3e-03 0.0e+00 1.0e-03 1.54e-07
## 16 20 -9.525e+02 2.78e-10 6.43e-11 1.9e-05 0.0e+00 1.0e-05 6.43e-11
## 17 21 -9.525e+02 -3.74e-12 1.07e-14 8.9e-07 0.0e+00 4.2e-07 1.07e-14
##
## ***** RELATIVE FUNCTION CONVERGENCE *****
##
## FUNCTION -9.524929e+02 RELDX 8.950e-07
## FUNC. EVALS 21 GRAD. EVALS 17
## PRELDF 1.074e-14 NPRELDF 1.074e-14
##
## I FINAL X(I) D(I) G(I)
##
## 1 4.938559e-02 1.000e+00 1.506e-04
## 2 2.207532e-01 1.000e+00 6.686e-05
## 3 1.849754e-01 1.000e+00 -1.711e-05fit$coe[3]## b1
## 0.1849754#so take the value of lambda as 0.83The two sets of lambda are shown below. The first list presents lambda values from 2018-2022, while the second presents values from 2013-2017.
lambda = c(0.83, 0.82, 0.81, 0.88, 0.74, 0.87, 0.84, 0.66, 0.63, 0.86, 0.79)
lambda2 = c(0.80, 0.71, 0.60, 0.59, 0.97, 0.17, 0.43, 0.98, 0.27, 0.97, 0.68)Section 4 - Ordinary least squares based algorithm for prediction
In this section I build a modified prediction algorithm (using ordinary least squares) to compare the Sharpe ratio of predicted daily returns of each splited dataset.
For the algorithm, it takes in a warmup time t0, a window length D, and the appropriately normalised 10+q return series as i.
I assume the process of validation is stable to use the lambda derived from Section 2 to form the further calculation instead of randomly selecting numbers from 0 to 1, which are respectively lambda = c(0.83, 0.82, 0.81, 0.88, 0.74, 0.87, 0.84, 0.66, 0.63, 0.86, 0.79) for period 2018-2022 and lambda2 = c(0.80, 0.71, 0.60, 0.59, 0.97, 0.17, 0.43, 0.98, 0.27, 0.97, 0.68) for period 2013-2017.
Explanations of functions
vol.exp.sm
The vol.exp.sm function performs exponential smoothing
on a given time series using a specified parameter lambda. The function
first calculates the square \(\sigma^2\) of the input time series data
x, then applies an exponential smoothing formula to
calculate the smoothed variance. This formula calculates the exponential
moving average of \(x^2\) and weights
it by \(1-\lambda\) to compute the
smoothed variance \(\sigma^2\).
The function then calculates the standard deviation \(\sigma\) by taking the square root of \(\sigma^2\). The residuals are computed by dividing the input data by the standard deviation \(\sigma\). If any of the residuals are NA, the function sets them to 0. The squared residuals are calculated as the square of the residuals.
The function returns a list containing the smoothed variance \(\sigma^2\), standard deviation \(\sigma\), residuals resid, and
squared residuals sq.resid.
first.acf.squares.train
The first.acf.squares.train function takes
x returned by the pred.footsie.prepare
function as input, and calculates the acfs (autocorrelation functions)
of the squared residuals after removing the volatility for each
covariate and the responses in the training set. This is done by
iterating through each element and computing the volatility using the
vol.exp.sm function with a specified parameter lambda.
The function then adds up the first acfs for each covariate and response to determine the amount of acf that has been removed. The goal is to choose lambda such that as much of the acf as possible is removed.
It also prepares the data by subtracting the volatility from each element in the training, validation, and testing sets. The squared residuals are then calculated and returned as a list along with the total sum of the first acfs, the modified training and validation data, and the modified testing data.
thresh.reg
The basic thresh.reg function performs a linear
regression of y onto the covariates in x that
have a marginal correlation with y greater than the
threshold th. The function returns the estimated
coefficients beta, the predicted value for a new
x if provided, and the condition number of the Gram matrix,
which is a measure of the stability of the regression.
The extended rolling.thresh.reg function performs
rolling window predictions over the training set using the exponential
smoothing parameter lambda and the thresh.reg
function. The function returns the mean squared error, predicted values,
true values, and condition number for each prediction.
The rolling.thresh.reg.valid and
rolling.thresh.reg.test function is the same as
rolling.thresh.reg, but for the validation set and the test
set respectively.
Define functions
Concept for the algorithm
For lag q = 0
Period 2018-2022
The algorithm was executed with a value of q = 0, and Sharpe ratios were computed for the training, validation, and test sets. The analysis was performed with D values ranging from 50 to 250, incremented by 20 each time.
The results are shown below.
#try for q=0
q=0
data = pred.footsie.prepare(q)
sc = sharpe.curves(data, lambda, 0, 250,
win = seq(from = 50, to = 250, by = 20))
sc$train.curve## [1] 0.288909429 0.005442172 0.126718369 1.143259894 1.818796966 1.282369271
## [7] 1.358573962 1.694284181 1.686846047 1.478086762 1.240343052sc$valid.curve## [1] 2.08515017 0.54744531 0.97678281 1.22540208 -1.38038619 0.05621205
## [7] -0.18647376 0.41126894 0.66305254 2.56760291 -0.73524306sc$test.curve## [1] 0.62740795 -0.82415943 -0.71262423 -1.97201833 0.03242541 -1.41874528
## [7] -4.36865332 -2.58540230 -0.89542739 -0.60868537 0.35550674Negative values were observed in the training and validation sets, indicating unoptimal performance on all three datasets. A rolling window approach was employed with a range of 90 to 250 and an interval of 16 (10% of the total range) to improve the Sharpe ratio of the training set. The starting value of 90 was identified as the range where all Sharpe ratios of the training set were positive. This approach aimed to enhance the strategy’s performance on the validation and test sets, and seek the chance to improve training set.
Finalised Sharpe ratio for period 1 when q=0
- Comparison between Sharpe ratios of two methods when q=0 (period 1): Back to the pcr method comparison
sc1 = sharpe.curves(data, lambda, 0, 250,
win = seq(from = 90, to = 250, by = 16))
sc1$train.curve## [1] 0.1267184 0.8772180 1.3138329 1.4789300 0.8265222 1.3585740 1.3659940
## [8] 1.6700413 1.7047476 1.8005283 1.2403431sc1$valid.curve## [1] 0.9767828 0.4099034 1.1266990 -2.1931209 1.9387690 -0.1864738
## [7] 0.4723983 -0.8592184 1.8303872 0.8934714 -0.7352431sc1$test.curve## [1] -0.7126242 -0.2620386 -2.0991942 -1.9946648 -1.1814791 -4.3686533
## [7] -1.8362397 -2.7116051 0.5856022 -0.8784990 0.3555067Although the Sharpe ratio for the training set improved, negative values were still observed in the validation set. Given the discrete nature of the dataset, the current window size was deemed to be sufficient, and reducing it further was considered unnecessary.
Period 2013-2017
Same algorithm was run with same q and starting window size for second period in order to make comparison.
q=0
data = pred.footsie.prepare2(q)
sc = sharpe.curves(data, lambda2, 0, 250,
win = seq(from = 50, to = 250, by = 20))
sc$train.curve## [1] 0.7269293 0.7202627 -0.2367161 -0.4703070 -0.7760157 -0.6834376
## [7] -0.6676578 -0.6644061 -0.3184242 -0.6634266 -1.1783759sc$valid.curve## [1] -0.37306456 -1.36767980 -1.65365232 -1.03702278 -0.90953472 0.39739909
## [7] 1.17494616 1.67507426 1.82184436 3.43778926 -0.03431468sc$test.curve## [1] -0.4976426 0.5784132 1.5141209 2.7637462 1.1698597 1.6626344
## [7] 1.0825274 0.7922051 0.6344409 -0.1342038 0.5777836Negative values were observed in nearly all of the training set and the first half of the validation set, while the test set performed relatively well. To increase the positive Sharpe ratios over the selected window, a new window with a range of 140 to 250 and an interval of 11 (10% of the total range) was chosen. The starting value of 140 was determined as the range where all Sharpe ratios of the validation set were positive. Given the poor performance of the training set, the focus was shifted to improving the performance of the other two sets.
Finalised Sharpe ratio for period 2 when q=0
- Comparison between Sharpe ratios of two methods when q=0 (period 2): Back to the pcr method comparison
sc2 = sharpe.curves(data, lambda2, 0, 250,
win = seq(from = 140, to = 250, by = 11))
sc2$train.curve## [1] -0.24506147 -0.33211342 -0.81554555 -1.07170372 -0.30874134 -0.84400289
## [7] 0.01072897 -0.86531656 -0.31378449 -1.19032519 -1.17837591sc2$valid.curve## [1] 0.34544957 -0.52987053 0.91569178 2.07014049 1.06568482 1.01920678
## [7] 0.58094377 1.74881228 0.64574701 1.18249995 -0.03431468sc2$test.curve## [1] 1.8696241 1.6626344 1.0966443 0.7095946 -0.5042961 1.1244827
## [7] 0.6344409 0.6344409 -0.4024697 -0.3704866 0.5777836Although the Sharpe ratio for the validation set improved, negative values were still observed all around the training set. Given the discrete nature of the dataset, the current window size was deemed to be sufficient, and reducing it further was considered unnecessary.
Concept for unsatisfactory Sharpe ratio
One of the reasons for having an unsatisfactory Sharpe ratio is multicollinearity. It means that covariates can be correlated with each other due to using too many covariates in linear regression. In the example above, this could be the case that they are all world indices, so their performances are usually highly correlated (see log-price plot in Section 2). In mathematical terms, multicollinearity means that Z has small eigenvalues.
The condition number of a matrix A is defined as \[ \text{Cond.number} = \frac{|\lambda|_{\text{max}}(A)}{|\lambda|_{\text{min}}(A)} \] where the denominator is the minimum absolute value of the eigenvalues of A, and the numerator is for the maximum. It could be an indicator of multicollinearity. The closer the number of covariates p is to the sample size, the larger the condition number is for a sample covariance matrix.
Plots of condition numbers when q=0
Period 2018-2022
par(mfrow=c(1,2))
plot.ts(sc1$condnum.test[4,])
plot.ts(sc1$condnum.test[11,])
The figure shows that the condition numbers are larger when D = 138, and
are better when D = 250. This is due to the dimension-to-sample size
ratio p/D becoming smaller. So overall, if D = 250 the performance could
be better.
Period 2013-2017
par(mfrow=c(1,2))
plot.ts(sc2$condnum.test[4,])
plot.ts(sc2$condnum.test[11,])
The figure shows that the condition numbers are larger when D = 173, and
are better when D = 250. This is due to the dimension-to-sample size
ratio p/D becoming smaller. So overall, if D = 250 the performance could
be better.
The condition number of two chosen windows from 2018 to 2022 is much smaller than the condition number from 2013 to 2017, which indicates that the factors affecting the system during 2018-2022 were less correlated compared to those during time 2013-2017.
For lag q = 1
Period 2018-2022
The algorithm was executed with a value of q = 1, and Sharpe ratios were computed for the training, validation, and test sets. The analysis was performed with D values ranging from 50 to 250, incremented by 20 each time.
The results are shown below.
#comparing for q=1
q=1
data = pred.footsie.prepare(q)
sc = sharpe.curves(data, lambda, 0, 250,
win = seq(from = 50, to = 250, by = 20))
sc$train.curve## [1] 0.05502265 0.70632808 0.47227334 1.02048456 1.92111410 1.47232150
## [7] 1.18914380 2.27279469 2.27543839 2.06599931 1.55933491sc$valid.curve## [1] -0.85831166 -0.56273112 2.33202645 0.09667942 0.53525205 1.85500677
## [7] 2.40466683 0.80954795 0.95614431 -1.01899719 -1.74796416sc$test.curve## [1] -2.03897188 -1.17671840 -0.04758179 -1.40999215 -0.53337453 -1.89383842
## [7] -5.14879249 -3.65680311 -0.70168638 -1.20057522 -0.76315020Initial analysis of the validation set revealed a lack of suitable ranges for the Sharpe ratio of the test set, as many negative values were observed when q=1 and the window was set to 50 to 250 with an interval of 20. However, the tail and end of the validation set were easier to discard.
Therefore, a new window with a range of 70 to 200 and an interval of 13 (1/10 of the total length) was chosen to obtain a better performing Sharpe ratio for both the training and validation sets, while disregarding the test set. Despite this improvement, some negative values were still observed.
Finalised Sharpe ratio for period 1 when q=1
- Comparison between Sharpe ratios of two methods when q=1 (period 1): Back to the pcr method comparison
sc3 = sharpe.curves(data, lambda, 0, 250,
win = seq(from = 70, to = 200, by = 13))
sc3$train.curve## [1] 0.7063281 0.4101290 0.5050730 1.0961117 1.7525484 2.1136680 1.5967829
## [8] 1.5629842 1.5554598 1.6384738 1.3908318sc3$valid.curve## [1] -0.5627311 -0.4825086 1.9367585 1.6308856 1.0863461 -0.3453824
## [7] 1.7201920 2.1082070 1.0634594 0.8095479 2.1176964sc3$test.curve## [1] -1.1767184 -1.2943346 -0.9502633 -1.6993592 -3.7637211 -0.6536850
## [7] -0.7065885 -2.2316506 -5.5857032 -2.8805391 -4.5783216Period 2013-2017
Same algorithm was run with same q and starting window size for second period in order to make comparison.
#comparing for q=1
q=1
data = pred.footsie.prepare2(q)
sc = sharpe.curves(data, lambda2, 0, 250,
win = seq(from = 50, to = 250, by = 20))
sc$train.curve## [1] 0.73956263 1.10876149 -0.01030753 -0.54581921 -0.16017773 0.48229287
## [7] 0.37586318 0.03900317 -0.05653700 -0.06619021 -0.22624015sc$valid.curve## [1] 8.071442e-01 -9.965492e-05 1.739491e+00 2.747935e+00 2.325690e+00
## [6] 1.747166e+00 1.564987e+00 1.396696e+00 -5.170489e-01 -3.635322e-02
## [11] 5.400142e-01sc$test.curve## [1] 1.3280169 -0.4012840 2.1631655 0.7506369 -2.3516580 -3.1276473
## [7] -2.7877433 -0.3189292 0.3660890 -0.5839908 -0.8673514Initial analysis of the training set revealed a lack of suitable ranges for the Sharpe ratio of the test set, as many negative values were observed when q=1 and the window was set to 50 to 250 with an interval of 20. However, the tail of the validation set were easier to discard. To further explore the suitable range of overall window, the Sharpe ratios of validation set indicates that the tail and end could be cut.
Therefore, a new window with a range of 80 to 200 and an interval of 12 (1/10 of the total length) was chosen to obtain a better performing Sharpe ratio for both the training and validation sets, while disregarding the test set.
Finalised Sharpe ratio for period 2 when q=1
- Comparison between Sharpe ratios of two methods when q=1 (period 2): Back to the pcr method comparison
sc4 = sharpe.curves(data, lambda2, 0, 250,
win = seq(from = 80, to = 200, by = 12))
sc4$train.curve## [1] -0.2259147 0.1073450 -1.1702863 -0.4055254 -0.4731658 -0.3240934
## [7] 0.1148436 1.0215571 0.2003118 -0.1433442 -0.1760901sc4$valid.curve## [1] 0.465236 1.922340 2.775324 2.364594 2.580534 1.925930 1.747166 1.257476
## [9] 1.678764 1.764121 0.675960sc4$test.curve## [1] -0.6205264 2.5327871 0.4292176 -1.3031001 -1.1563914 -2.5112523
## [7] -3.0330409 -2.3305672 -2.0196486 -0.3189292 -0.8441028Plots of condition numbers when q=1
Period 2018-2022
Then plot the condition numbers of different D to find which value of D is better.
par(mfrow=c(1,2))
plot.ts(sc3$condnum.test[4,])
plot.ts(sc3$condnum.test[11,])
The difference between condition numbers is quite significant between D
= 109 and D = 200, so if D = 200 the performance could be better.
Period 2013-2017
par(mfrow=c(1,2))
plot.ts(sc4$condnum.test[4,])
plot.ts(sc4$condnum.test[11,])
The difference between condition numbers is quite significant between D
= 116 and D = 200, so if D = 200 the performance could be better.
Short summary of each lag’s performance
In summary, the algorithm’s performance was greatly improved by adjusting the window for each value, particularly when q=1 in both period. Despite the better performance of the q=0 setting in terms of multicollinearity (as indicated by a smaller conditional number at D=250), the appropriateness of the algorithm was ultimately determined by the Sharpe ratio of the validation set.
Section 5 - Principal component regression based algorithm for prediction
In order to improve upon ordinary least squares, the one-day-ahead S&P500 return is to be predicted using the factors from the 10 stocks chosen as covariates using technique called principal component regression.
Instead of determining the number of factors using a scree plot for each window, it treats the number of factors as another tuning parameter, on top of the window length. And in each window, perform a multi-factor analysis, use the estimated factor series as the covariates and a linear model for predicting the one-day-ahead S&500 return.
Then output Sharpe ratios should then be dependent on window length as well as number of factors considered in each window.
Period 2018-2022
For lag q = 0
The object data.0.dev is a new data that removes the
effect of volatility.
#when q=0
q=0
data = pred.footsie.prepare(q)
data.0.dev = first.acf.squares.train(data, lambda)And then use the function prcomp() in basic R to find the principal component analysis of each set. Then summary them to read the proportion of variance.
#For varying volatility removed data
train.pcr = prcomp(data.0.dev$x.train.dev)
valid.pcr = prcomp(data.0.dev$x.valid.dev)
test.pcr = prcomp(data.0.dev$x.test.dev)
summary(train.pcr)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 2.6842 1.3664 1.3168 1.19752 1.08130 1.02419 0.85255
## Proportion of Variance 0.4379 0.1135 0.1054 0.08716 0.07106 0.06375 0.04418
## Cumulative Proportion 0.4379 0.5514 0.6568 0.74392 0.81498 0.87874 0.92291
## PC8 PC9 PC10
## Standard deviation 0.7876 0.64670 0.47951
## Proportion of Variance 0.0377 0.02542 0.01397
## Cumulative Proportion 0.9606 0.98603 1.00000summary(valid.pcr)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 2.6396 1.1759 1.11519 0.98771 0.92514 0.84880 0.75740
## Proportion of Variance 0.5079 0.1008 0.09066 0.07111 0.06239 0.05252 0.04182
## Cumulative Proportion 0.5079 0.6087 0.69932 0.77044 0.83283 0.88534 0.92716
## PC8 PC9 PC10
## Standard deviation 0.67137 0.58615 0.45271
## Proportion of Variance 0.03286 0.02504 0.01494
## Cumulative Proportion 0.96002 0.98506 1.00000summary(test.pcr)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 3.0014 1.2695 1.14461 0.96459 0.84869 0.72012 0.69178
## Proportion of Variance 0.5885 0.1053 0.08559 0.06078 0.04706 0.03388 0.03126
## Cumulative Proportion 0.5885 0.6938 0.77939 0.84017 0.88723 0.92111 0.95237
## PC8 PC9 PC10
## Standard deviation 0.56852 0.46726 0.43306
## Proportion of Variance 0.02112 0.01426 0.01225
## Cumulative Proportion 0.97348 0.98775 1.00000The first two principal component contains the majority of the data set fitted. (more than 50%) So it is acceptable if the model only consider number of factors up to 2.
Define function
The algorithm is quite similar to the one in Section 4, but I still
need to change three rolling.thresh.reg functions. As I only want look
up to 2 factors, so I create a new definition PCR_analysis
which only contains the first and second factor. The function would only
be run when PCR_analysis is FALSE.
Comparison between Sharpe ratios of two methods when q=0 (period 1)
Then use the same adjusted window for q=0 in Section 4 to see if there is any improvements.
sc5 = sharpe.curves(data, lambda, 0, 250,
win = seq(from = 90, to = 250, by = 16))
sc5$train.curve## [1] 0.1757264 0.9757615 1.5664493 2.0600731 1.6556380 2.0875776 2.3280796
## [8] 1.6126724 1.9027072 1.8989932 1.5689709sc5$valid.curve## [1] -0.7586326 -0.3827531 -0.6138769 -1.4611453 0.5585362 -0.1212978
## [7] 0.4723983 -0.7979058 0.9394031 0.8934714 -1.3672700sc5$test.curve## [1] -0.51567264 -1.43506828 -3.31979610 -1.99477101 -1.18148847 -3.40529525
## [7] -1.83612699 -2.71151378 0.58563135 -1.16307247 0.07267749As the model is intended for real-world prediction, the test set is of greater value in assessing the accuracy of the predicted data.
The performance of the validation set is superior to that of the ordinary least squares method. However, the most different Sharpe ratios in the end of the test set suggest that the investment is likely to generate more negative returns for the same level of risk, whereas the other ratios indicate little difference between the two methods.
For lag q = 1
Then use the same adjusted window for q=1 in Section 4 to see if there is any improvements.
#when q=1
q=1
data = pred.footsie.prepare(q)
data.1.dev = first.acf.squares.train(data, lambda)
#For varying volatility removed data
train.pcr = prcomp(data.1.dev$x.train.dev)
valid.pcr = prcomp(data.1.dev$x.valid.dev)
test.pcr = prcomp(data.1.dev$x.test.dev)Comparison between Sharpe ratios of two methods when q=1 (period 1)
sc6 = sharpe.curves(data, lambda, 0, 250,
win = seq(from = 70, to = 200, by = 13))
sc6$train.curve## [1] 0.9356807 0.5640615 0.5310563 0.8152106 1.6440209 1.9164978 1.5324106
## [8] 1.6344186 1.4140016 1.6657542 1.5381076sc6$valid.curve## [1] -0.5627311 -0.4825086 1.9367585 1.6308856 1.0863461 -0.3453824
## [7] 1.7201920 2.1082070 1.0634594 0.8095479 2.1176964sc6$test.curve## [1] -1.0520645 -2.7256508 -1.4341069 -1.4184541 -3.7637201 -0.5659552
## [7] -2.0238918 -1.1989365 -4.7949615 -3.3289347 -4.5287404The performance of test set is significantly worse than using original method. Hence it could be concluded that using principle component regression would not improve the correctness and liability of the prediction.
Period 2013-2017
Same approach is applied again for the interval 2013-2017, so I would only focus on results for this part.
For lag q = 0
#when q=0
q=0
data = pred.footsie.prepare2(q)
data.0.dev.1 = first.acf.squares.train(data, lambda2)#For varying volatility removed data
train.pcr= prcomp(data.0.dev.1$x.train.dev)
valid.pcr = prcomp(data.0.dev.1$x.valid.dev)
test.pcr= prcomp(data.0.dev.1$x.test.dev)
summary(train.pcr)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 2.8554 2.0472 1.7451 1.50446 1.37788 1.11014 1.0520
## Proportion of Variance 0.3404 0.1749 0.1271 0.09449 0.07926 0.05145 0.0462
## Cumulative Proportion 0.3404 0.5153 0.6424 0.73693 0.81619 0.86763 0.9138
## PC8 PC9 PC10
## Standard deviation 0.9826 0.86579 0.59075
## Proportion of Variance 0.0403 0.03129 0.01457
## Cumulative Proportion 0.9541 0.98543 1.00000summary(valid.pcr)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 3.0775 2.0672 1.7113 1.2672 1.14762 1.03783 0.92663
## Proportion of Variance 0.4147 0.1871 0.1282 0.0703 0.05766 0.04716 0.03759
## Cumulative Proportion 0.4147 0.6018 0.7300 0.8003 0.85797 0.90513 0.94272
## PC8 PC9 PC10
## Standard deviation 0.80272 0.6157 0.53355
## Proportion of Variance 0.02821 0.0166 0.01246
## Cumulative Proportion 0.97094 0.9875 1.00000summary(test.pcr)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 3.2954 2.5787 2.2974 1.48882 1.23587 1.01832 0.94239
## Proportion of Variance 0.3628 0.2221 0.1763 0.07405 0.05102 0.03464 0.02967
## Cumulative Proportion 0.3628 0.5849 0.7612 0.83528 0.88630 0.92094 0.95061
## PC8 PC9 PC10
## Standard deviation 0.81991 0.69022 0.57434
## Proportion of Variance 0.02246 0.01591 0.01102
## Cumulative Proportion 0.97307 0.98898 1.00000The first two principal component contains the majority of the data set fitted. (more than 50%) So it is acceptable if the model only consider number of factors up to 2.
Comparison between Sharpe ratios of two methods when q=0 (period 2)
sc7 = sharpe.curves(data, lambda2, 0, 250,
win = seq(from = 140, to = 250, by = 11))
sc7$train.curve## [1] -0.2819389 -0.5025671 -0.9816810 -0.8535206 -0.4089415 -0.6268541
## [7] -0.1987614 -0.5749441 -0.8264982 -1.0042689 -1.4691360sc7$valid.curve## [1] 0.8597323 -0.1331840 1.3151252 2.5270877 1.0252911 1.4437199
## [7] 0.9788364 2.1544802 3.7657200 1.5835556 0.4789857sc7$test.curve## [1] 1.6354878 1.6485152 1.0966446 -0.8361755 0.4302447 1.5605065
## [7] 0.3030262 -0.1341960 -0.1341960 -0.3704802 0.3412968The training set’s overall performance is inferior, but the mean Sharpe ratio of the test sets is higher on average.
For lag q = 1
#when q=1
q=1
data = pred.footsie.prepare2(q)
data.1.dev.1 = first.acf.squares.train(data, lambda2)#For varying volatility removed data
train.pcr= prcomp(data.1.dev.1$x.train.dev)
valid.pcr = prcomp(data.1.dev.1$x.valid.dev)
test.pcr= prcomp(data.1.dev.1$x.test.dev)
summary(train.pcr)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 3.4489 2.5961 1.77123 1.46408 1.30743 1.16932 1.11859
## Proportion of Variance 0.3788 0.2147 0.09992 0.06827 0.05444 0.04355 0.03985
## Cumulative Proportion 0.3788 0.5935 0.69342 0.76169 0.81613 0.85968 0.89953
## PC8 PC9 PC10 PC11
## Standard deviation 1.01752 0.94697 0.89895 0.6437
## Proportion of Variance 0.03298 0.02856 0.02574 0.0132
## Cumulative Proportion 0.93250 0.96106 0.98680 1.0000summary(valid.pcr)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 3.3034 2.0763 1.9348 1.55707 1.28524 1.11772 0.9972
## Proportion of Variance 0.3918 0.1548 0.1344 0.08705 0.05931 0.04485 0.0357
## Cumulative Proportion 0.3918 0.5466 0.6810 0.76801 0.82732 0.87217 0.9079
## PC8 PC9 PC10 PC11
## Standard deviation 0.8815 0.85363 0.79488 0.65461
## Proportion of Variance 0.0279 0.02616 0.02268 0.01539
## Cumulative Proportion 0.9358 0.96193 0.98461 1.00000summary(test.pcr)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 3.589 2.8272 2.1034 1.46815 1.25902 1.12479 1.0067
## Proportion of Variance 0.380 0.2358 0.1305 0.06359 0.04676 0.03732 0.0299
## Cumulative Proportion 0.380 0.6158 0.7463 0.80994 0.85670 0.89403 0.9239
## PC8 PC9 PC10 PC11
## Standard deviation 0.93267 0.80060 0.74052 0.72067
## Proportion of Variance 0.02566 0.01891 0.01618 0.01532
## Cumulative Proportion 0.94959 0.96850 0.98468 1.00000Comparison between Sharpe ratios of two methods when q=1 (period 2)
sc8 = sharpe.curves(data, lambda2, 0, 250,
win = seq(from = 140, to = 250, by = 11))
sc8$train.curve## [1] 0.08511143 0.63104257 0.82157329 0.36851057 0.17690784 0.29786746
## [7] -0.12719183 -0.50871600 -0.08605264 -0.41360598 -0.09340578sc8$valid.curve## [1] 2.5876839 2.2462173 0.8201173 2.1737285 1.9993262 1.2633651 1.2166221
## [8] 2.1697203 1.4271187 2.3361972 0.3420592sc8$test.curve## [1] -4.0080949 -3.1505013 -2.0518056 -2.9627858 -1.0077529 -0.6561661
## [7] -0.8235817 -0.2096081 -0.5839908 0.1324463 -0.8673514The Sharpe ratios in the tail of the test set suggest that the investment is likely to generate more negative returns for the same level of risk.
Short summary of each lag’s performance
From the above comparison, it is evident that after implementing the PCR method, the Sharpe ratios for all periods significantly decreased, and more negative values appeared. Therefore, a simple conclusion can be drawn that the OLS-based algorithm is more appropriate for the model analysis presented in this report.
Section 6 - Study for Risk measures
The model constructed earlier was found to be less accurate, prompting the development of a risk management model for further analysis and comparison of data from both periods.
For the risk management analysis, the most well-known measure - Value-at-Risk (VaR) is used.
Value-at-risk (VaR) with all its associated problems remains the workhorse model in risk management. One key problem of VaR is that it does not properly account for volatility clustering, which means that VaR limits are breached in serial dependence across time. As a result, risk is underestimated during a crisis. A powerful approach to solving this problem is to combine VaR with GARCH models, which take conditional volatility into account. To illustrate this method, a GARCH(1,1) model is applied.
Importing data - 1
To generate log-price plots for each stock, the data used in the
previous step were transformed into text files. However, in order to
generate Chartseries plots for each stock, it was necessary to repeat
the importation process to obtain xts objects.
Using ARIMA model to estimate Value-at-Risk
The following section only focus on the method of deriving VaR, so here only contain an example of S&p500.
rets = diff(GSPC[,4]) / GSPC[,4][-length(GSPC[,4])]
rets <- na.omit(rets)
model.arima = auto.arima(rets , max.order = c(3 , 0 ,3) , stationary = TRUE , trace = T , ic = 'aicc')##
## Fitting models using approximations to speed things up...
##
## ARIMA(2,0,2) with non-zero mean : -7100
## ARIMA(0,0,0) with non-zero mean : -7052.698
## ARIMA(1,0,0) with non-zero mean : -7091.263
## ARIMA(0,0,1) with non-zero mean : -7084.178
## ARIMA(0,0,0) with zero mean : -7054.412
## ARIMA(1,0,2) with non-zero mean : -7102.506
## ARIMA(0,0,2) with non-zero mean : -7104.92
## ARIMA(0,0,3) with non-zero mean : -7103.516
## ARIMA(1,0,1) with non-zero mean : -7096.347
## ARIMA(1,0,3) with non-zero mean : -7101.393
## ARIMA(0,0,2) with zero mean : -7106.617
## ARIMA(0,0,1) with zero mean : -7085.775
## ARIMA(1,0,2) with zero mean : -7104.207
## ARIMA(0,0,3) with zero mean : -7105.199
## ARIMA(1,0,1) with zero mean : -7097.973
## ARIMA(1,0,3) with zero mean : -7103.109
##
## Now re-fitting the best model(s) without approximations...
##
## ARIMA(0,0,2) with zero mean : -7106.562
##
## Best model: ARIMA(0,0,2) with zero meanThe chosen model for the return is ARIMA(0,0,2) with a zero mean, selected based on the AIC scores computed for various ARIMA models. This model represents a 2nd order Moving Average (MA(2)) model.
Model estimation - 1
To estimate the parameter coefficients, Maximum Likelihood method was employed. For the selected ARIMA(0, 0, 2) model, the estimated coefficients are as follows:
model.arima## Series: rets
## ARIMA(0,0,2) with zero mean
##
## Coefficients:
## ma1 ma2
## -0.1614 0.1395
## s.e. 0.0280 0.0286
##
## sigma^2 = 0.0001858: log likelihood = 3556.29
## AIC=-7106.58 AICc=-7106.56 BIC=-7091.22Therefore the process can be described as:
\[ r_t = -0.1614 \times r_{t-1}+0.1395 \times r_{t-2} + \epsilon_t \]
Using GARCH(1,1) to estimate Value-at-Risk
Model estimation - 2
In order to test the Heteroscedasticity of the squared residuals, significance testing on a1 and β1 parameters is applied.
ar.res = model.arima$residuals
model.spec = ugarchspec(variance.model = list(model = 'sGARCH' , garchOrder = c(1 , 1)) ,
mean.model = list(armaOrder = c(0 , 0)))
model.fit = ugarchfit(spec = model.spec , data = ar.res , solver = 'solnp')
options(scipen = 999)
model.fit@fit$matcoef## Estimate Std. Error t value Pr(>|t|)
## mu 0.000966497563 0.000237837552 4.063688 0.00004830342069728211
## omega 0.000005448641 0.000003118223 1.747355 0.08057582566843146665
## alpha1 0.254625652941 0.033861825831 7.519549 0.00000000000005506706
## beta1 0.736417539351 0.033423555877 22.032890 0.00000000000000000000Both a1 and β1 are significantly different from zero, therefore it is reasonable to assume time-varying volatility of the residuals.
he GARCH equation can be expressed as: \[ \sigma^2= 0.0000054+0.2546\epsilon^2_{t−1}+0.7364\sigma^2_{t−1} \] Given that 0 < β1 < 1, as lag increases the effect of the squared residual decreases.
Deriving VaR
The 95% VaR could be estimated below, where red bars refer to returns lower than 5% quantile.
quantile(rets, 0.05)## 5%
## -0.02164038qplot(rets, geom = 'histogram', bins = 30) +
geom_histogram(aes(x = rets[rets < quantile(rets, 0.05)]), fill = 'red', bins = 30) +
labs(x = 'Daily Returns')
In order to model more adequately the thickness of tails, other distributional assumptions could be used for stock returns. The t-distribution is symmetric and bell-shaped, like the normal distribution, but has heavier tails, meaning that it is more prone to producing values that fall far from its mean. We use the fitdist function from rugarch package to get the fitting parameters of t-distribution.
fitdist(distribution = 'std' , x = rets)$pars## mu sigma shape
## 0.001026709 0.016484844 2.544354797cat("For a = 0.05 the quantile value of normal distribution is: " ,
qnorm(p = 0.05) , "\n" ,
"For a = 0.05 the quantile value of t-distribution is: " ,
qdist(distribution = 'std' , shape = 2.544354797 , p = 0.05) , "\n" , "\n" ,
'For a = 0.01 the quantile value of normal distribution is: ' ,
qnorm(p = 0.01) , "\n" ,
"For a = 0.01 the quantile value of t-distribution is: " ,
qdist(distribution = 'std' , shape = 2.544354797 , p = 0.01) , sep = "")## For a = 0.05 the quantile value of normal distribution is: -1.644854
## For a = 0.05 the quantile value of t-distribution is: -1.172791
##
## For a = 0.01 the quantile value of normal distribution is: -2.326348
## For a = 0.01 the quantile value of t-distribution is: -2.432742qplot(y = rets , x = 1:1236 , geom = 'point') + geom_point(colour = 'lightgrey' , size = 2) +
geom_line(aes(y = model.fit@fit$sigma*(-1.172791) , x = 1:1236) , colour = 'red') +
labs(x = '' , y = 'Daily Returns' , title = 'Value at Risk Visualisation')+ theme_light() +
theme(legend.position = 'none') 
Red line denotes VaR produced by GARCH model
VaR forecasting
The ugarchroll method allows to perform a rolling estimation and forecasting of a model combination. It returns the distributional forecast parameters necessary to calculate any required measure on the forecasted density. I set the 75% of the data as test set and perform a rolling moving 1-step ahead forecast of the conditional standard deviation, \(\hat σ_{t+1|t}\), and re-estimate GARCH parameters every 50 observations.
model.roll = ugarchroll(spec = model.spec , data =rets , n.start = 309 , refit.every = 50 , refit.window = 'moving')
VaR95_td = mean(rets) + model.roll@forecast$density[,'Sigma']*qdist(distribution='std', shape=2.544354797, p=0.05)
qplot(y = VaR95_td , x = 1:927 , geom = 'line') +
geom_point(aes(x = 1:927 , y = rets[310:1236] , color = as.factor(rets[310:1236] < VaR95_td)) , size = 1.5) + scale_color_manual(values = c('gray' , 'red')) +
labs(y = 'Daily Returns' , x = 'Test set Observation') + theme_light() +
theme(legend.position = 'none')
cat('\n' , 'Number of exceptions with GARCH approach: ' , (sum(rets < VaR95_td)) , sep = '')##
## Number of exceptions with GARCH approach: 153Black line represent the daily forecasted VaR given by the GARCH model and red points refer to returns lower than VaR, and the last step is to count the number of exceptions and compare it with the one generated with delta-normal approach.
Section 7 - Conclusion
The study has shown that the analysis of the daily closing prices of selected stocks and the S&P 500 index can assist potential investors or risk managers in making investment decisions. The use of statistical and financial theorems, including ordinary least square analysis and principal component regression analysis, to construct the model, and exponential smoothing to improve accuracy in capturing volatility, has proven effective in predicting returns, Value at Risk (VaR).
Reference
How to use R prcomp results for prediction? https://stats.stackexchange.com/questions/72839/how-to-use-r-prcomp-results-for-prediction
Exponential smoothing versus GARCH (1,1) for conditional variance https://stats.stackexchange.com/questions/229685/exponential-smoothing-versus-garch1-1-for-conditional-variance
Principal Components Analysis (PCA) in R https://www.benjaminbell.co.uk/2018/02/principal-components-analysis-pca-in-r.html
Value at Risk (VaR) in R Programming Language https://maung-sutikno.medium.com/value-at-risk-var-in-r-programming-language-2c28cc92fbc2
Value at risk estimation using GARCH model https://rpubs.com/ionaskel/VaR_Garch_market_risk
VaR with GARCH(1,1) http://rstudio-pubs-static.s3.amazonaws.com/241528_8cb5ef038484439cb106f3ef883329ac.html
Appendix
Text files-1
- Export data into text files-1: Back to the main body
#export the data to text files
GSPC = as.data.frame(GSPC)
GSPC = cbind(as.numeric(as.Date(rownames(GSPC))), GSPC)
write.table(GSPC, "GSPC.txt", row.names=FALSE)
AAPL = as.data.frame(AAPL)
AAPL = cbind(as.numeric(as.Date(rownames(AAPL))), AAPL)
write.table(AAPL, "AAPL.txt", row.names=FALSE)
AMZN = as.data.frame(AMZN)
AMZN = cbind(as.numeric(as.Date(rownames(AMZN))), AMZN)
write.table(AMZN, "AMZN.txt", row.names=FALSE)
AVGO = as.data.frame(AVGO)
AVGO = cbind(as.numeric(as.Date(rownames(AVGO))), AVGO)
write.table(AVGO, "AVGO.txt", row.names=FALSE)
COST = as.data.frame(COST)
COST = cbind(as.numeric(as.Date(rownames(COST))), COST)
write.table(COST, "COST.txt", row.names=FALSE)
GOOG = as.data.frame(GOOG)
GOOG = cbind(as.numeric(as.Date(rownames(GOOG))), GOOG)
write.table(GOOG, "GOOG.txt", row.names=FALSE)
META = as.data.frame(META)
META = cbind(as.numeric(as.Date(rownames(META))), META)
write.table(META, "META.txt", row.names=FALSE)
MSFT = as.data.frame(MSFT)
MSFT = cbind(as.numeric(as.Date(rownames(MSFT))), MSFT)
write.table(MSFT, "MSFT.txt", row.names=FALSE)
NVDA = as.data.frame(NVDA)
NVDA = cbind(as.numeric(as.Date(rownames(NVDA))), NVDA)
write.table(NVDA, "NVDA.txt", row.names=FALSE)
PEP = as.data.frame(PEP)
PEP = cbind(as.numeric(as.Date(rownames(PEP))), PEP)
write.table(PEP, "PEP.txt", row.names=FALSE)
TSLA = as.data.frame(TSLA)
TSLA = cbind(as.numeric(as.Date(rownames(TSLA))), TSLA)
write.table(TSLA, "TSLA.txt", row.names=FALSE)Text files-2
- Export data into text files-2: Back to the main body
#export the data to text files
GSPC2 = as.data.frame(GSPC2)
GSPC2 = cbind(as.numeric(as.Date(rownames(GSPC2))), GSPC2)
write.table(GSPC2, "GSPC2.txt", row.names=FALSE)
AAPL2 = as.data.frame(AAPL2)
AAPL2 = cbind(as.numeric(as.Date(rownames(AAPL2))), AAPL2)
write.table(AAPL2, "AAPL2.txt", row.names=FALSE)
AMZN2 = as.data.frame(AMZN2)
AMZN2 = cbind(as.numeric(as.Date(rownames(AMZN2))), AMZN2)
write.table(AMZN2, "AMZN2.txt", row.names=FALSE)
AVGO2 = as.data.frame(AVGO2)
AVGO2 = cbind(as.numeric(as.Date(rownames(AVGO2))), AVGO2)
write.table(AVGO2, "AVGO2.txt", row.names=FALSE)
COST2 = as.data.frame(COST2)
COST2 = cbind(as.numeric(as.Date(rownames(COST2))), COST2)
write.table(COST2, "COST2.txt", row.names=FALSE)
GOOG2 = as.data.frame(GOOG2)
GOOG2 = cbind(as.numeric(as.Date(rownames(GOOG2))), GOOG2)
write.table(GOOG2, "GOOG2.txt", row.names=FALSE)
META2 = as.data.frame(META2)
META2 = cbind(as.numeric(as.Date(rownames(META2))), META2)
write.table(META2, "META2.txt", row.names=FALSE)
MSFT2 = as.data.frame(MSFT2)
MSFT2 = cbind(as.numeric(as.Date(rownames(MSFT2))), MSFT2)
write.table(MSFT2, "MSFT2.txt", row.names=FALSE)
NVDA2 = as.data.frame(NVDA2)
NVDA2 = cbind(as.numeric(as.Date(rownames(NVDA2))), NVDA2)
write.table(NVDA2, "NVDA2.txt", row.names=FALSE)
PEP2 = as.data.frame(PEP2)
PEP2 = cbind(as.numeric(as.Date(rownames(PEP2))), PEP2)
write.table(PEP2, "PEP2.txt", row.names=FALSE)
TSLA2 = as.data.frame(TSLA2)
TSLA2 = cbind(as.numeric(as.Date(rownames(TSLA2))), TSLA2)
write.table(TSLA2, "TSLA2.txt", row.names=FALSE)Function of log-price plot of 2018-2022
- Log-price plot of 2018-2022: Back to the main body
#define read.bossa.data function
read.bossa.data <- function(vec.names) {
p <- length(vec.names)
n1 <- 20000
dates <- matrix(99999999, p, n1)
closes <- matrix(0, p, n1)
max.n2 <- 0
for (i in 1:p) {
filename <- paste("C:/Users/liugu/Desktop/ST326 Project/", vec.names[i], ".txt", sep="")
tmp <- scan(filename, list(date=numeric(), NULL, NULL, NULL,
close=numeric(), NULL, NULL), skip=1, sep="")
n2 <- length(tmp$date)
max.n2 <- max(n2, max.n2)
dates[i,1:n2] <- tmp$date
closes[i,1:n2] <- tmp$close
}
dates <- dates[,1:max.n2]
closes <- closes[,1:max.n2]
days <- rep(0, n1)
arranged.closes <- matrix(0, p, n1)
date.indices <- starting.indices <- rep(1, p)
already.started <- rep(0, p)
day <- 1
while(max(date.indices) <= max.n2) {
current.dates <- current.closes <- rep(0, p)
for (i in 1:p) {
current.dates[i] <- dates[i,date.indices[i]]
current.closes[i] <- closes[i,date.indices[i]]
}
min.indices <- which(current.dates == min(current.dates))
days[day] <- current.dates[min.indices[1]]
arranged.closes[min.indices,day] <- log(current.closes[min.indices])
arranged.closes[-min.indices,day] <- arranged.closes[-min.indices, max(day-1, 1)]
already.started[min.indices] <- 1
starting.indices[-which(already.started == 1)] <- starting.indices[-which(already.started == 1)] + 1
day <- day + 1
date.indices[min.indices] <- date.indices[min.indices] + 1
}
days <- days[1:(day-1)]
arranged.closes <- arranged.closes[,1:(day-1)]
max.st.ind <- max(starting.indices)
r <- matrix(0, p, (day-max.st.ind-1))
for (i in 1:p) {
r[i,] <- diff(arranged.closes[i,max.st.ind:(day-1)])
r[i,] <- r[i,] / sqrt(var(r[i,]))
r[i,r[i,]==0] <- rnorm(sum(r[i,]==0))
}
return(list(dates=dates, closes=closes, days=days, arranged.closes=arranged.closes, starting.indices=starting.indices, r=r))
}Define function of volatility
- Function of volatility: Back to the main body
#since the company selected all have the closing values on the same day, so set all shift indices as 0.(looking at data up t--the current time)
pred.footsie.prepare <- function(max.lag = 5, split = c(50, 25), mask = rep(1, 10)) {
ind <- read.bossa.data(c("AAPL","MSFT","AMZN","GOOG","AVGO","TSLA","NVDA","PEP","META","COST","GSPC"))
d <- dim(ind$r)
start.index <- max(3, max.lag + 1)
y <- matrix(0, d[2] - start.index + 1, 1)
x <- matrix(0, d[2] - start.index + 1, d[1] - 1 + max.lag)
y[,1] <- ind$r[1,start.index:d[2]]
for (i in 1:max.lag) {
x[,i] <- ind$r[1,(start.index-i):(d[2]-i)]
}
shift.indices <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
for (i in 2:(d[1])) {
x[,i+max.lag-1] <- ind$r[i,(start.index-1-shift.indices[i-1]):(d[2]-1-shift.indices[i-1])]
}
end.training <- round(split[1] / 100 * d[2])
end.validation <- round(sum(split[1:2]) / 100 * d[2])
x <- x[,as.logical(c(rep(1, max.lag), mask))]
y.train <- as.matrix(y[1:end.training], end.training, 1)
x.train <- x[1:end.training,]
y.valid <- as.matrix(y[(end.training+1):(end.validation)], end.validation-end.training, 1)
x.valid <- x[(end.training+1):(end.validation),]
y.test <- as.matrix(y[(end.validation+1):(d[2] - start.index + 1)], d[2]-start.index-end.validation+1, 1)
x.test <- x[(end.validation+1):(d[2] - start.index + 1),]
list(x=x, y=y, x.train=x.train, y.train=y.train, x.valid=x.valid, y.valid=y.valid, x.test=x.test, y.test=y.test)
}
pred.footsie.prepare2 <- function(max.lag = 5, split = c(50, 25), mask = rep(1, 10)) {
ind <- read.bossa.data(c("AAPL2","MSFT2","AMZN2","GOOG2","AVGO2","TSLA2","NVDA2","PEP2","META2","COST2","GSPC2"))
d <- dim(ind$r)
start.index <- max(3, max.lag + 1)
y <- matrix(0, d[2] - start.index + 1, 1)
x <- matrix(0, d[2] - start.index + 1, d[1] - 1 + max.lag)
y[,1] <- ind$r[1,start.index:d[2]]
for (i in 1:max.lag) {
x[,i] <- ind$r[1,(start.index-i):(d[2]-i)]
}
shift.indices <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
for (i in 2:(d[1])) {
x[,i+max.lag-1] <- ind$r[i,(start.index-1-shift.indices[i-1]):(d[2]-1-shift.indices[i-1])]
}
end.training <- round(split[1] / 100 * d[2])
end.validation <- round(sum(split[1:2]) / 100 * d[2])
x <- x[,as.logical(c(rep(1, max.lag), mask))]
y.train <- as.matrix(y[1:end.training], end.training, 1)
x.train <- x[1:end.training,]
y.valid <- as.matrix(y[(end.training+1):(end.validation)], end.validation-end.training, 1)
x.valid <- x[(end.training+1):(end.validation),]
y.test <- as.matrix(y[(end.validation+1):(d[2] - start.index + 1)], d[2]-start.index-end.validation+1, 1)
x.test <- x[(end.validation+1):(d[2] - start.index + 1),]
list(x=x, y=y, x.train=x.train, y.train=y.train, x.valid=x.valid, y.valid=y.valid, x.test=x.test, y.test=y.test)
}Detailed explanation of concept for using exponential smoothing to estimate daily volatilities
- Concept for using exponential smoothing to estimate daily volatilities: Back to the main body
When considering the simple exponential smoothing model: \[ \sigma_t^2 = \lambda \sigma^2_{t-1}+(1-\lambda)\ x^2_{t-1} \] And GARCH (1,1) in its simplest form (w.r.t. the conditional mean equation): \[ r_t = \sigma_t \varepsilon_t;\\ \varepsilon_t \sim i.i.d(0,1);\\ \sigma^2_t = a_0+b_0x^2_{t-1}+c_0\sigma^2_{t-1} \] And rewrite the distribution so that it becomes: \[ (1-c_0B)\sigma^2_t = a_0+b_0x^2_{t-1},\text{so}\\ \sigma^2_t = (1-c_0B)^{-1}(a_0+b_0x^2_{t-1})\\ = \sum_{k\ge0}(c_0B)^{k}(a_0+b_0x^2_{t-1})\\ = \frac{a_0}{1-c_0}+b_0\sum_{k\ge0}c_0^{k}x^2_{t-k-1} \] The mean of the volatility is defined by: \[ E(\sigma^2_t) = \frac{a_0}{1-b_0-c_0} \] And this mean is large since 1-b0-c0 ≈ 0, so the constant a0/(1−c0) is tiny compared to the mean. The constant term could be ignored and only consider the rest elements. Moreover, \[ b_0(1+c_0+c_0^2+...) = \frac{b_0}{1-c_0} \approx \frac{1-c_0}{1-c_0} = 1 \] Hence, we can approximately model the volatility using an exponential filter with 0 < λ < 1. (Take b0=1-λ and c0=λ.) \[ \sigma^2_t = (1-\lambda)\sum_{k\ge0}\lambda^kx^2_{t-k-1} \] And this could be simplified into the form which is similar to the simple exponential smoothing model: \[ \sigma^2_t = (1-\lambda)x^2_{t-1}+\lambda\sigma^2_{t-1} \] Since it is easy to show that \[ E(\sigma^2_t) = E(\sigma^2_{t-1}) \] And \[ \tilde{\sigma}^2_{t+k}= \sigma^2_t \] can be used as a forecast when current time is t and k > 0 is not too large.
Detailed code of two period lambda
- Define functions: Back to the main body
#λ for period 2018-2022
#Apple Inc
lp.AAPL = log(AAPL[,5])
r1 = diff(lp.AAPL)
r1 = r1[2:length(r1)]
fit = garch(r1)##
## ***** ESTIMATION WITH ANALYTICAL GRADIENT *****
##
##
## I INITIAL X(I) D(I)
##
## 1 4.026443e-04 1.000e+00
## 2 5.000000e-02 1.000e+00
## 3 5.000000e-02 1.000e+00
##
## IT NF F RELDF PRELDF RELDX STPPAR D*STEP NPRELDF
## 0 1 -4.173e+03
## 1 7 -4.174e+03 2.94e-04 4.68e-04 2.1e-04 4.3e+09 2.1e-05 1.00e+06
## 2 8 -4.174e+03 1.55e-05 1.77e-05 2.0e-04 2.0e+00 2.1e-05 1.12e+01
## 3 15 -4.186e+03 2.91e-03 4.63e-03 3.9e-01 2.0e+00 6.3e-02 1.11e+01
## 4 18 -4.211e+03 5.88e-03 5.20e-03 6.9e-01 1.9e+00 2.5e-01 9.39e-01
## 5 20 -4.217e+03 1.40e-03 1.36e-03 7.5e-02 2.0e+00 5.1e-02 1.39e+02
## 6 22 -4.230e+03 3.10e-03 2.90e-03 1.2e-01 2.0e+00 1.0e-01 1.52e+04
## 7 23 -4.250e+03 4.84e-03 6.58e-03 1.8e-01 2.0e+00 2.0e-01 1.95e+06
## 8 35 -4.259e+03 1.99e-03 4.67e-03 1.5e-05 3.1e+00 2.0e-05 9.23e-01
## 9 36 -4.259e+03 1.05e-04 7.82e-05 1.5e-05 2.0e+00 2.0e-05 4.82e-01
## 10 37 -4.259e+03 1.35e-05 1.60e-05 1.5e-05 2.0e+00 2.0e-05 6.46e-01
## 11 44 -4.264e+03 1.15e-03 2.36e-03 5.8e-02 2.0e+00 8.3e-02 6.30e-01
## 12 45 -4.266e+03 4.49e-04 4.88e-04 1.5e-02 0.0e+00 2.4e-02 4.88e-04
## 13 47 -4.269e+03 6.50e-04 3.86e-04 2.5e-02 0.0e+00 4.7e-02 3.86e-04
## 14 49 -4.270e+03 2.99e-04 3.37e-04 1.7e-02 1.2e+00 3.4e-02 9.57e-04
## 15 50 -4.271e+03 1.13e-04 1.52e-04 1.7e-02 5.4e-01 3.4e-02 1.87e-04
## 16 51 -4.271e+03 1.19e-06 4.82e-06 2.7e-03 0.0e+00 5.2e-03 4.82e-06
## 17 52 -4.271e+03 1.55e-06 2.11e-06 1.2e-03 0.0e+00 2.0e-03 2.11e-06
## 18 54 -4.271e+03 8.47e-10 1.53e-07 8.2e-05 1.4e+00 1.4e-04 1.72e-07
## 19 55 -4.271e+03 1.11e-07 4.31e-08 3.1e-05 0.0e+00 5.8e-05 4.31e-08
## 20 66 -4.271e+03 -7.03e-15 1.00e-16 2.8e-15 7.4e+04 4.7e-15 1.92e-10
##
## ***** FALSE CONVERGENCE *****
##
## FUNCTION -4.270735e+03 RELDX 2.809e-15
## FUNC. EVALS 66 GRAD. EVALS 20
## PRELDF 1.004e-16 NPRELDF 1.923e-10
##
## I FINAL X(I) D(I) G(I)
##
## 1 1.920489e-05 1.000e+00 9.160e+01
## 2 1.264059e-01 1.000e+00 6.444e-03
## 3 8.333217e-01 1.000e+00 -2.215e-02fit$coe[3]## b1
## 0.8333217#so take the value of lambda as 0.83
#Amazon.com Inc
lp.AMZN = log(AMZN[,5])
r2 = diff(lp.AMZN)
r2 = r2[2:length(r2)]
fit = garch(r2)##
## ***** ESTIMATION WITH ANALYTICAL GRADIENT *****
##
##
## I INITIAL X(I) D(I)
##
## 1 4.591286e-04 1.000e+00
## 2 5.000000e-02 1.000e+00
## 3 5.000000e-02 1.000e+00
##
## IT NF F RELDF PRELDF RELDX STPPAR D*STEP NPRELDF
## 0 1 -4.087e+03
## 1 7 -4.089e+03 3.05e-04 5.02e-04 2.6e-04 3.1e+09 2.6e-05 7.69e+05
## 2 8 -4.089e+03 1.13e-05 1.27e-05 2.5e-04 2.0e+00 2.6e-05 1.86e+01
## 3 15 -4.114e+03 6.25e-03 1.00e-02 5.1e-01 2.0e+00 1.1e-01 1.84e+01
## 4 18 -4.161e+03 1.12e-02 8.34e-03 7.9e-01 1.9e+00 4.2e-01 7.37e-01
## 5 20 -4.168e+03 1.79e-03 2.19e-03 8.1e-02 2.0e+00 8.5e-02 4.50e+02
## 6 21 -4.191e+03 5.29e-03 8.18e-03 7.0e-02 2.0e+00 8.5e-02 2.14e+02
## 7 22 -4.200e+03 2.13e-03 3.03e-03 6.1e-02 2.0e+00 8.5e-02 4.18e+01
## 8 30 -4.200e+03 1.25e-04 2.68e-04 2.5e-06 4.7e+00 3.6e-06 1.31e-01
## 9 31 -4.200e+03 3.26e-06 3.08e-06 2.4e-06 2.0e+00 3.6e-06 2.15e-01
## 10 32 -4.200e+03 3.96e-08 6.20e-08 2.4e-06 2.0e+00 3.6e-06 1.97e-01
## 11 38 -4.200e+03 1.43e-05 1.94e-05 2.5e-03 2.0e+00 3.7e-03 1.97e-01
## 12 41 -4.201e+03 3.03e-04 1.67e-04 1.5e-02 0.0e+00 2.7e-02 1.67e-04
## 13 42 -4.203e+03 3.59e-04 6.53e-04 5.3e-02 3.3e-01 1.1e-01 7.03e-04
## 14 43 -4.203e+03 9.28e-05 2.44e-04 1.5e-02 0.0e+00 2.6e-02 2.44e-04
## 15 44 -4.204e+03 4.52e-05 3.80e-05 4.2e-03 0.0e+00 8.2e-03 3.80e-05
## 16 45 -4.204e+03 3.45e-06 5.68e-06 4.0e-03 1.6e-01 8.2e-03 5.78e-06
## 17 46 -4.204e+03 3.32e-07 4.70e-07 7.0e-04 0.0e+00 1.2e-03 4.70e-07
## 18 47 -4.204e+03 6.49e-08 2.55e-08 2.1e-04 0.0e+00 3.8e-04 2.55e-08
## 19 63 -4.204e+03 -2.60e-15 3.59e-15 8.7e-15 3.9e+05 1.4e-14 2.19e-09
##
## ***** FALSE CONVERGENCE *****
##
## FUNCTION -4.203533e+03 RELDX 8.702e-15
## FUNC. EVALS 63 GRAD. EVALS 19
## PRELDF 3.595e-15 NPRELDF 2.191e-09
##
## I FINAL X(I) D(I) G(I)
##
## 1 2.001557e-05 1.000e+00 1.074e+03
## 2 1.610626e-01 1.000e+00 1.722e-01
## 3 8.087633e-01 1.000e+00 1.880e-01fit$coe[3]## b1
## 0.8087633#so take the value of lambda as 0.81
#Microsoft Corp
lp.MSFT = log(MSFT[,5])
r3 = diff(lp.MSFT)
r3 = r3[2:length(r3)]
fit = garch(r3)##
## ***** ESTIMATION WITH ANALYTICAL GRADIENT *****
##
##
## I INITIAL X(I) D(I)
##
## 1 3.470341e-04 1.000e+00
## 2 5.000000e-02 1.000e+00
## 3 5.000000e-02 1.000e+00
##
## IT NF F RELDF PRELDF RELDX STPPAR D*STEP NPRELDF
## 0 1 -4.294e+03
## 1 7 -4.299e+03 1.32e-03 2.12e-03 4.2e-04 5.2e+09 4.2e-05 5.51e+06
## 2 8 -4.299e+03 3.19e-05 4.50e-05 4.1e-04 2.0e+00 4.2e-05 4.29e+01
## 3 9 -4.299e+03 9.40e-06 9.55e-06 4.2e-04 2.0e+00 4.2e-05 4.24e+01
## 4 15 -4.318e+03 4.37e-03 6.01e-03 3.0e-01 2.0e+00 4.3e-02 4.21e+01
## 5 18 -4.334e+03 3.69e-03 3.53e-03 6.2e-01 2.0e+00 1.7e-01 9.09e-01
## 6 20 -4.380e+03 1.03e-02 7.76e-03 4.3e-01 2.0e+00 3.4e-01 8.76e+01
## 7 22 -4.390e+03 2.46e-03 2.55e-03 5.7e-02 2.0e+00 6.9e-02 2.08e+04
## 8 24 -4.417e+03 6.05e-03 6.00e-03 9.7e-02 2.0e+00 1.4e-01 7.21e+02
## 9 31 -4.419e+03 3.84e-04 9.12e-04 2.7e-06 1.8e+01 4.1e-06 1.43e+02
## 10 32 -4.419e+03 6.17e-06 5.55e-06 2.7e-06 2.0e+00 4.1e-06 3.32e+01
## 11 33 -4.419e+03 2.69e-07 3.47e-07 2.7e-06 2.0e+00 4.1e-06 3.43e+01
## 12 41 -4.422e+03 6.94e-04 2.47e-03 4.2e-02 2.0e+00 6.8e-02 3.42e+01
## 13 42 -4.423e+03 2.12e-04 1.09e-03 2.1e-02 0.0e+00 3.5e-02 1.09e-03
## 14 43 -4.424e+03 1.85e-04 1.25e-04 2.3e-03 0.0e+00 4.0e-03 1.25e-04
## 15 44 -4.424e+03 4.40e-05 4.95e-05 2.4e-03 1.2e+00 4.0e-03 9.17e-05
## 16 45 -4.424e+03 1.29e-05 1.54e-05 2.4e-03 4.5e-01 4.0e-03 1.72e-05
## 17 46 -4.424e+03 2.17e-07 4.22e-07 3.5e-04 0.0e+00 6.1e-04 4.22e-07
## 18 47 -4.424e+03 5.06e-08 1.94e-08 8.5e-05 0.0e+00 1.9e-04 1.94e-08
## 19 48 -4.424e+03 9.26e-08 4.19e-08 3.0e-04 0.0e+00 6.6e-04 4.19e-08
## 20 49 -4.424e+03 7.59e-09 4.57e-09 8.7e-05 0.0e+00 1.9e-04 4.57e-09
## 21 60 -4.424e+03 -3.70e-15 1.26e-15 3.6e-15 3.3e+05 5.9e-15 2.50e-10
##
## ***** FALSE CONVERGENCE *****
##
## FUNCTION -4.423863e+03 RELDX 3.626e-15
## FUNC. EVALS 60 GRAD. EVALS 21
## PRELDF 1.259e-15 NPRELDF 2.505e-10
##
## I FINAL X(I) D(I) G(I)
##
## 1 1.177111e-05 1.000e+00 9.406e+02
## 2 1.597299e-01 1.000e+00 8.374e-02
## 3 8.165350e-01 1.000e+00 1.683e-01fit$coe[3]## b1
## 0.816535#so take the value of lambda as 0.82
#Alphabet Inc
lp.GOOG = log(GOOG[,5])
r4 = diff(lp.GOOG)
r4 = r4[2:length(r4)]
fit = garch(r4)##
## ***** ESTIMATION WITH ANALYTICAL GRADIENT *****
##
##
## I INITIAL X(I) D(I)
##
## 1 3.526665e-04 1.000e+00
## 2 5.000000e-02 1.000e+00
## 3 5.000000e-02 1.000e+00
##
## IT NF F RELDF PRELDF RELDX STPPAR D*STEP NPRELDF
## 0 1 -4.249e+03
## 1 7 -4.250e+03 2.01e-04 3.07e-04 1.5e-04 6.0e+09 1.5e-05 9.27e+05
## 2 8 -4.250e+03 1.72e-05 1.93e-05 1.3e-04 2.0e+00 1.5e-05 9.79e+00
## 3 15 -4.260e+03 2.19e-03 2.95e-03 3.0e-01 2.0e+00 4.4e-02 9.73e+00
## 4 18 -4.268e+03 1.92e-03 1.76e-03 6.2e-01 1.9e+00 1.7e-01 2.17e-01
## 5 20 -4.293e+03 5.77e-03 4.21e-03 4.3e-01 2.0e+00 3.5e-01 3.46e+01
## 6 22 -4.299e+03 1.57e-03 1.54e-03 5.7e-02 2.0e+00 7.0e-02 1.05e+04
## 7 23 -4.304e+03 1.21e-03 3.15e-03 9.7e-02 2.0e+00 1.4e-01 3.40e+05
## 8 31 -4.306e+03 2.84e-04 2.21e-03 6.7e-06 3.4e+00 1.0e-05 8.95e+00
## 9 32 -4.307e+03 3.83e-04 3.02e-04 6.2e-06 2.0e+00 1.0e-05 5.78e+00
## 10 33 -4.307e+03 1.83e-05 2.43e-05 6.6e-06 2.0e+00 1.0e-05 1.02e+00
## 11 34 -4.307e+03 1.15e-06 9.53e-07 6.7e-06 2.0e+00 1.0e-05 1.53e+00
## 12 40 -4.308e+03 3.77e-05 6.59e-05 6.9e-03 2.0e+00 1.1e-02 1.51e+00
## 13 42 -4.309e+03 3.19e-04 1.68e-04 1.6e-02 0.0e+00 3.3e-02 1.68e-04
## 14 44 -4.311e+03 3.94e-04 4.41e-04 2.7e-02 1.3e+00 5.3e-02 1.64e-03
## 15 45 -4.311e+03 8.06e-05 3.11e-04 2.6e-02 9.1e-01 5.3e-02 5.90e-04
## 16 54 -4.311e+03 3.28e-05 9.11e-05 6.8e-07 2.8e+00 1.2e-06 1.20e-04
## 17 55 -4.311e+03 1.04e-05 1.05e-05 5.1e-07 2.0e+00 1.2e-06 7.93e-05
## 18 56 -4.311e+03 9.90e-09 3.93e-08 5.2e-07 2.0e+00 1.2e-06 6.08e-05
## 19 64 -4.311e+03 2.29e-05 3.78e-05 4.8e-03 7.7e-01 1.1e-02 6.05e-05
## 20 67 -4.311e+03 1.70e-08 2.11e-08 4.0e-05 1.8e+00 7.7e-05 7.31e-08
## 21 69 -4.311e+03 1.05e-09 4.03e-09 1.6e-05 1.8e+00 3.0e-05 3.75e-08
## 22 72 -4.311e+03 4.12e-10 3.33e-10 1.6e-06 2.0e+00 3.1e-06 3.27e-08
## 23 83 -4.311e+03 -4.85e-15 2.34e-18 2.3e-15 1.2e+03 4.1e-15 3.34e-08
##
## ***** FALSE CONVERGENCE *****
##
## FUNCTION -4.311317e+03 RELDX 2.281e-15
## FUNC. EVALS 83 GRAD. EVALS 23
## PRELDF 2.338e-18 NPRELDF 3.345e-08
##
## I FINAL X(I) D(I) G(I)
##
## 1 1.715719e-05 1.000e+00 2.407e+00
## 2 7.571116e-02 1.000e+00 4.131e-01
## 3 8.803744e-01 1.000e+00 -2.790e-01fit$coe[3]## b1
## 0.8803744#so take the value of lambda as 0.88
#Tesla Inc
lp.TSLA = log(TSLA[,5])
r5 = diff(lp.TSLA)
r5 = r5[2:length(r5)]
fit = garch(r5)##
## ***** ESTIMATION WITH ANALYTICAL GRADIENT *****
##
##
## I INITIAL X(I) D(I)
##
## 1 1.516962e-03 1.000e+00
## 2 5.000000e-02 1.000e+00
## 3 5.000000e-02 1.000e+00
##
## IT NF F RELDF PRELDF RELDX STPPAR D*STEP NPRELDF
## 0 1 -3.338e+03
## 1 7 -3.338e+03 6.25e-05 1.19e-04 4.0e-04 2.5e+08 4.0e-05 1.50e+04
## 2 13 -3.342e+03 1.12e-03 1.67e-03 2.8e-01 2.0e+00 4.1e-02 2.84e+00
## 3 16 -3.353e+03 3.26e-03 2.88e-03 5.7e-01 2.0e+00 1.6e-01 6.45e-01
## 4 17 -3.371e+03 5.26e-03 7.23e-03 4.2e-01 2.0e+00 3.3e-01 7.08e+01
## 5 25 -3.374e+03 9.23e-04 4.63e-03 1.3e-04 2.8e+00 1.4e-04 2.05e-02
## 6 26 -3.376e+03 6.52e-04 4.89e-04 1.2e-04 2.0e+00 1.4e-04 2.43e-03
## 7 27 -3.376e+03 5.10e-05 6.44e-05 1.3e-04 2.0e+00 1.4e-04 6.93e-03
## 8 28 -3.376e+03 4.37e-06 4.48e-06 1.3e-04 2.0e+00 1.4e-04 5.75e-03
## 9 34 -3.378e+03 4.69e-04 1.11e-03 8.2e-02 1.5e+00 9.7e-02 5.80e-03
## 10 35 -3.380e+03 7.46e-04 4.22e-04 2.2e-02 0.0e+00 3.3e-02 4.22e-04
## 11 36 -3.386e+03 1.61e-03 1.80e-03 7.5e-02 8.1e-01 1.3e-01 3.12e-03
## 12 37 -3.389e+03 7.55e-04 1.05e-03 7.6e-02 6.4e-01 1.3e-01 1.41e-03
## 13 39 -3.389e+03 1.13e-04 2.31e-04 4.8e-03 1.5e+00 8.5e-03 1.19e-03
## 14 40 -3.389e+03 6.87e-05 4.00e-04 4.1e-03 9.5e-01 8.5e-03 7.00e-04
## 15 41 -3.390e+03 1.15e-04 1.35e-04 4.4e-03 1.7e+00 8.5e-03 5.22e-04
## 16 42 -3.390e+03 3.20e-05 3.86e-05 4.5e-03 9.2e-01 8.5e-03 8.28e-05
## 17 43 -3.390e+03 2.13e-05 1.57e-05 3.0e-03 0.0e+00 5.4e-03 1.57e-05
## 18 44 -3.390e+03 1.70e-06 8.15e-07 8.7e-04 0.0e+00 1.6e-03 8.15e-07
## 19 45 -3.390e+03 1.63e-07 2.44e-08 1.7e-04 0.0e+00 3.1e-04 2.44e-08
## 20 46 -3.390e+03 5.73e-09 2.56e-10 8.2e-06 0.0e+00 1.9e-05 2.56e-10
## 21 47 -3.390e+03 -9.82e-10 1.23e-11 8.7e-07 0.0e+00 2.1e-06 1.23e-11
##
## ***** RELATIVE FUNCTION CONVERGENCE *****
##
## FUNCTION -3.389706e+03 RELDX 8.679e-07
## FUNC. EVALS 47 GRAD. EVALS 21
## PRELDF 1.230e-11 NPRELDF 1.230e-11
##
## I FINAL X(I) D(I) G(I)
##
## 1 9.593162e-05 1.000e+00 -1.864e+01
## 2 7.576528e-02 1.000e+00 -4.609e-02
## 3 8.668286e-01 1.000e+00 -4.011e-02fit$coe[3]## b1
## 0.8668286#so take the value of lambda as 0.87
#NVIDIA Corp
lp.NVDA = log(NVDA[,5])
r6 = diff(lp.NVDA)
r6 = r6[2:length(r6)]
fit = garch(r6)##
## ***** ESTIMATION WITH ANALYTICAL GRADIENT *****
##
##
## I INITIAL X(I) D(I)
##
## 1 9.681559e-04 1.000e+00
## 2 5.000000e-02 1.000e+00
## 3 5.000000e-02 1.000e+00
##
## IT NF F RELDF PRELDF RELDX STPPAR D*STEP NPRELDF
## 0 1 -3.625e+03
## 1 6 -3.625e+03 3.28e-05 8.66e-04 1.0e-03 3.1e+08 1.0e-04 1.36e+05
## 2 7 -3.626e+03 2.06e-04 2.38e-04 5.0e-04 2.0e+00 5.0e-05 1.02e+01
## 3 14 -3.640e+03 3.82e-03 6.42e-03 4.6e-01 2.0e+00 8.6e-02 1.03e+01
## 4 17 -3.655e+03 4.23e-03 3.67e-03 6.8e-01 1.9e+00 2.4e-01 3.99e-01
## 5 19 -3.697e+03 1.13e-02 9.92e-03 4.5e-01 2.0e+00 4.8e-01 1.36e+02
## 6 24 -3.708e+03 2.90e-03 5.75e-03 3.1e-05 2.8e+00 4.8e-05 6.35e-01
## 7 25 -3.709e+03 4.14e-04 1.36e-03 2.8e-05 2.0e+00 4.8e-05 2.65e+00
## 8 26 -3.710e+03 2.83e-04 2.20e-04 3.1e-05 2.0e+00 4.8e-05 1.08e+00
## 9 27 -3.710e+03 1.79e-05 2.27e-05 3.1e-05 2.0e+00 4.8e-05 6.14e-01
## 10 33 -3.711e+03 2.08e-04 3.97e-04 1.1e-02 2.0e+00 1.8e-02 6.76e-01
## 11 35 -3.713e+03 6.19e-04 3.58e-04 2.1e-02 0.0e+00 4.1e-02 3.58e-04
## 12 43 -3.713e+03 5.13e-06 9.85e-06 7.3e-07 8.1e+00 1.1e-06 1.61e-03
## 13 44 -3.713e+03 2.75e-09 2.95e-08 5.7e-07 2.0e+00 1.1e-06 1.22e-03
## 14 54 -3.715e+03 3.56e-04 4.25e-04 1.9e-02 7.5e-01 3.7e-02 1.22e-03
## 15 55 -3.715e+03 1.53e-04 1.75e-04 1.7e-02 3.7e-01 3.7e-02 1.92e-04
## 16 56 -3.715e+03 6.12e-06 1.65e-05 3.9e-03 0.0e+00 6.8e-03 1.65e-05
## 17 57 -3.715e+03 3.92e-06 3.91e-06 7.2e-04 0.0e+00 1.3e-03 3.91e-06
## 18 58 -3.715e+03 6.28e-08 2.04e-08 5.9e-05 0.0e+00 1.2e-04 2.04e-08
## 19 71 -3.715e+03 -4.77e-15 1.85e-15 6.0e-15 1.3e+06 1.0e-14 2.71e-09
##
## ***** FALSE CONVERGENCE *****
##
## FUNCTION -3.715397e+03 RELDX 6.043e-15
## FUNC. EVALS 71 GRAD. EVALS 19
## PRELDF 1.846e-15 NPRELDF 2.711e-09
##
## I FINAL X(I) D(I) G(I)
##
## 1 3.727150e-05 1.000e+00 6.751e+02
## 2 1.304261e-01 1.000e+00 8.947e-02
## 3 8.407792e-01 1.000e+00 2.206e-01fit$coe[3]## b1
## 0.8407792#so take the value of lambda as 0.84
#PepsiCo Inc
lp.PEP = log(PEP[,5])
r7 = diff(lp.PEP)
r7 = r7[2:length(r7)]
fit = garch(r7)##
## ***** ESTIMATION WITH ANALYTICAL GRADIENT *****
##
##
## I INITIAL X(I) D(I)
##
## 1 1.800762e-04 1.000e+00
## 2 5.000000e-02 1.000e+00
## 3 5.000000e-02 1.000e+00
##
## IT NF F RELDF PRELDF RELDX STPPAR D*STEP NPRELDF
## 0 1 -4.746e+03
## 1 7 -4.762e+03 3.26e-03 4.11e-03 2.7e-04 2.7e+10 2.7e-05 5.56e+07
## 2 8 -4.764e+03 4.97e-04 1.25e-03 2.7e-04 2.2e+00 2.7e-05 1.23e+02
## 3 9 -4.766e+03 2.68e-04 3.67e-04 2.4e-04 2.0e+00 2.7e-05 1.10e+02
## 4 10 -4.766e+03 2.08e-05 1.90e-05 2.6e-04 2.0e+00 2.7e-05 1.12e+02
## 5 17 -4.818e+03 1.09e-02 1.76e-02 4.5e-01 2.0e+00 8.2e-02 1.11e+02
## 6 20 -4.867e+03 9.96e-03 1.20e-02 5.7e-01 2.0e+00 1.7e-01 1.08e+01
## 7 27 -4.868e+03 2.05e-04 4.50e-04 1.2e-05 1.8e+02 5.5e-06 1.11e-01
## 8 28 -4.868e+03 8.93e-07 8.84e-07 1.2e-05 2.0e+00 5.5e-06 3.80e-02
## 9 37 -4.906e+03 7.81e-03 1.49e-02 4.3e-01 9.7e-01 3.6e-01 3.78e-02
## 10 50 -4.918e+03 2.41e-03 9.09e-03 1.2e-05 2.6e+00 1.4e-05 1.20e-02
## 11 62 -4.924e+03 1.22e-03 7.79e-04 6.3e-02 0.0e+00 8.1e-02 7.79e-04
## 12 70 -4.926e+03 5.09e-04 7.23e-04 2.2e-06 4.3e+00 2.9e-06 1.90e-03
## 13 71 -4.926e+03 8.41e-06 1.19e-05 2.1e-06 2.0e+00 2.9e-06 1.14e-04
## 14 72 -4.926e+03 7.21e-07 5.94e-07 2.1e-06 2.0e+00 2.9e-06 8.10e-05
## 15 73 -4.926e+03 8.65e-09 1.31e-08 2.1e-06 2.0e+00 2.9e-06 8.10e-05
## 16 80 -4.927e+03 5.33e-05 6.69e-05 1.6e-02 5.2e-01 2.2e-02 8.10e-05
## 17 81 -4.927e+03 9.02e-05 4.82e-05 9.7e-03 0.0e+00 1.6e-02 4.82e-05
## 18 83 -4.929e+03 3.39e-04 3.58e-04 3.4e-02 2.4e-01 6.5e-02 6.20e-04
## 19 84 -4.929e+03 8.49e-06 1.47e-04 3.0e-02 1.3e-01 6.5e-02 1.48e-04
## 20 85 -4.929e+03 6.94e-05 1.03e-04 1.5e-02 1.7e-01 3.3e-02 1.05e-04
## 21 86 -4.929e+03 8.83e-07 5.29e-06 6.4e-04 0.0e+00 1.3e-03 5.29e-06
## 22 87 -4.929e+03 2.38e-06 2.24e-06 6.4e-04 1.4e+00 1.3e-03 2.35e-06
## 23 88 -4.929e+03 5.51e-08 1.66e-07 8.2e-04 0.0e+00 1.7e-03 1.66e-07
## 24 89 -4.929e+03 1.31e-08 8.62e-10 3.7e-05 0.0e+00 5.9e-05 8.62e-10
## 25 90 -4.929e+03 -3.28e-09 4.42e-11 9.5e-06 0.0e+00 1.5e-05 4.42e-11
##
## ***** RELATIVE FUNCTION CONVERGENCE *****
##
## FUNCTION -4.929099e+03 RELDX 9.471e-06
## FUNC. EVALS 90 GRAD. EVALS 25
## PRELDF 4.425e-11 NPRELDF 4.425e-11
##
## I FINAL X(I) D(I) G(I)
##
## 1 1.160405e-05 1.000e+00 4.071e+01
## 2 1.337195e-01 1.000e+00 -3.271e-02
## 3 7.861464e-01 1.000e+00 -2.455e-02fit$coe[3]## b1
## 0.7861464#so take the value of lambda as 0.66
#Meta Platforms Inc
lp.META = log(META[,5])
r8 = diff(lp.META)
r8 = r8[2:length(r8)]
fit = garch(r8)##
## ***** ESTIMATION WITH ANALYTICAL GRADIENT *****
##
##
## I INITIAL X(I) D(I)
##
## 1 7.026494e-04 1.000e+00
## 2 5.000000e-02 1.000e+00
## 3 5.000000e-02 1.000e+00
##
## IT NF F RELDF PRELDF RELDX STPPAR D*STEP NPRELDF
## 0 1 -3.817e+03
## 1 7 -3.817e+03 9.89e-05 1.73e-04 2.3e-04 1.3e+09 2.3e-05 1.10e+05
## 2 8 -3.817e+03 2.94e-06 3.08e-06 2.2e-04 2.0e+00 2.3e-05 8.74e+00
## 3 15 -3.833e+03 4.18e-03 6.28e-03 4.8e-01 2.0e+00 9.3e-02 8.69e+00
## 4 18 -3.854e+03 5.55e-03 6.09e-03 7.5e-01 1.9e+00 3.7e-01 5.40e-01
## 5 20 -3.863e+03 2.17e-03 2.10e-03 7.9e-02 2.0e+00 7.3e-02 6.90e-01
## 6 22 -3.868e+03 1.45e-03 1.51e-03 6.8e-02 2.0e+00 7.3e-02 1.22e+00
## 7 28 -3.869e+03 7.46e-06 1.44e-05 2.2e-06 2.9e+01 2.5e-06 6.62e-03
## 8 29 -3.869e+03 5.19e-08 1.02e-07 2.2e-06 2.0e+00 2.5e-06 5.88e-03
## 9 37 -3.873e+03 1.06e-03 1.36e-03 6.6e-02 1.4e+00 8.2e-02 5.87e-03
## 10 38 -3.873e+03 1.45e-04 2.22e-04 2.3e-02 0.0e+00 3.6e-02 2.22e-04
## 11 39 -3.874e+03 1.35e-04 1.04e-04 2.3e-02 0.0e+00 3.3e-02 1.04e-04
## 12 40 -3.874e+03 6.66e-05 6.24e-05 1.9e-02 1.9e-02 3.3e-02 6.24e-05
## 13 41 -3.874e+03 1.71e-05 1.40e-05 9.2e-03 0.0e+00 1.3e-02 1.40e-05
## 14 42 -3.874e+03 2.04e-06 2.09e-06 2.8e-03 0.0e+00 3.6e-03 2.09e-06
## 15 44 -3.874e+03 6.27e-09 4.45e-08 1.8e-04 8.1e-01 2.5e-04 6.64e-08
## 16 64 -3.874e+03 -2.82e-15 2.44e-15 1.4e-14 2.3e+06 1.8e-14 2.24e-08
##
## ***** FALSE CONVERGENCE *****
##
## FUNCTION -3.874014e+03 RELDX 1.400e-14
## FUNC. EVALS 64 GRAD. EVALS 16
## PRELDF 2.439e-15 NPRELDF 2.238e-08
##
## I FINAL X(I) D(I) G(I)
##
## 1 1.239365e-04 1.000e+00 5.395e+02
## 2 2.845577e-01 1.000e+00 -1.614e-01
## 3 6.252760e-01 1.000e+00 -8.693e-03fit$coe[3]## b1
## 0.625276#so take the value of lambda as 0.63
#Costco Wholesale Corp
lp.COST = log(COST[,5])
r9 = diff(lp.COST)
r9 = r9[2:length(r9)]
fit = garch(r9)##
## ***** ESTIMATION WITH ANALYTICAL GRADIENT *****
##
##
## I INITIAL X(I) D(I)
##
## 1 2.092837e-04 1.000e+00
## 2 5.000000e-02 1.000e+00
## 3 5.000000e-02 1.000e+00
##
## IT NF F RELDF PRELDF RELDX STPPAR D*STEP NPRELDF
## 0 1 -4.575e+03
## 1 7 -4.576e+03 2.53e-04 3.80e-04 1.0e-04 1.7e+10 1.0e-05 3.31e+06
## 2 8 -4.576e+03 2.30e-05 2.65e-05 8.9e-05 2.0e+00 1.0e-05 1.35e+01
## 3 9 -4.576e+03 1.25e-06 1.26e-06 9.9e-05 2.0e+00 1.0e-05 1.35e+01
## 4 16 -4.587e+03 2.35e-03 3.13e-03 2.9e-01 2.0e+00 4.1e-02 1.34e+01
## 5 20 -4.641e+03 1.16e-02 7.49e-03 8.6e-01 1.8e+00 6.6e-01 3.11e-01
## 6 32 -4.665e+03 5.04e-03 1.65e-02 1.2e-05 3.4e+00 1.7e-05 2.69e-02
## 7 42 -4.677e+03 2.62e-03 2.23e-03 3.1e-02 1.8e+00 4.5e-02 9.67e-02
## 8 48 -4.680e+03 6.45e-04 1.01e-03 2.5e-06 1.1e+01 3.7e-06 1.33e-01
## 9 49 -4.680e+03 1.00e-05 1.42e-05 2.4e-06 2.0e+00 3.7e-06 1.89e-01
## 10 50 -4.680e+03 1.13e-06 1.06e-06 2.5e-06 2.0e+00 3.7e-06 1.82e-01
## 11 59 -4.692e+03 2.53e-03 3.28e-03 3.9e-02 1.9e+00 6.1e-02 1.81e-01
## 12 61 -4.694e+03 3.56e-04 4.67e-04 1.5e-02 9.0e-01 2.5e-02 7.35e-04
## 13 62 -4.694e+03 4.36e-05 7.49e-05 5.5e-03 0.0e+00 1.3e-02 7.49e-05
## 14 63 -4.694e+03 6.37e-05 6.02e-05 5.8e-03 3.4e-01 1.3e-02 6.09e-05
## 15 64 -4.694e+03 8.78e-06 1.30e-05 1.9e-03 0.0e+00 3.4e-03 1.30e-05
## 16 65 -4.694e+03 1.30e-06 1.38e-06 4.1e-04 0.0e+00 7.7e-04 1.38e-06
## 17 85 -4.694e+03 -2.52e-15 1.20e-14 9.5e-15 1.5e+05 1.6e-14 1.03e-09
##
## ***** FALSE CONVERGENCE *****
##
## FUNCTION -4.694159e+03 RELDX 9.464e-15
## FUNC. EVALS 85 GRAD. EVALS 17
## PRELDF 1.199e-14 NPRELDF 1.029e-09
##
## I FINAL X(I) D(I) G(I)
##
## 1 7.795296e-06 1.000e+00 3.472e+03
## 2 1.146403e-01 1.000e+00 3.305e-01
## 3 8.560145e-01 1.000e+00 4.610e-01fit$coe[3]## b1
## 0.8560145#so take the value of lambda as 0.86
#Broadcom Inc
lp.AVGO = log(AVGO[,5])
r10 = diff(lp.AVGO)
r10 = r10[2:length(r10)]
fit = garch(r10)##
## ***** ESTIMATION WITH ANALYTICAL GRADIENT *****
##
##
## I INITIAL X(I) D(I)
##
## 1 5.107930e-04 1.000e+00
## 2 5.000000e-02 1.000e+00
## 3 5.000000e-02 1.000e+00
##
## IT NF F RELDF PRELDF RELDX STPPAR D*STEP NPRELDF
## 0 1 -4.047e+03
## 1 6 -4.049e+03 6.77e-04 3.27e-03 1.0e-03 1.3e+09 1.0e-04 2.16e+06
## 2 7 -4.051e+03 4.46e-04 6.38e-04 8.8e-04 2.0e+00 1.0e-04 4.71e+01
## 3 8 -4.051e+03 4.59e-05 4.17e-05 9.8e-04 2.0e+00 1.0e-04 4.83e+01
## 4 14 -4.084e+03 7.96e-03 1.32e-02 4.6e-01 2.0e+00 8.6e-02 4.80e+01
## 5 15 -4.096e+03 2.93e-03 3.72e-03 3.2e-01 2.0e+00 8.6e-02 2.35e+00
## 6 17 -4.115e+03 4.68e-03 6.83e-03 3.9e-01 2.0e+00 2.0e-01 2.35e+00
## 7 18 -4.132e+03 3.98e-03 7.79e-03 2.2e-01 1.9e+00 2.0e-01 2.14e-01
## 8 19 -4.137e+03 1.30e-03 6.03e-03 1.5e-01 1.9e+00 2.0e-01 7.27e-02
## 9 20 -4.139e+03 3.65e-04 1.70e-02 1.5e-01 1.3e+00 2.0e-01 1.71e-02
## 10 21 -4.139e+03 1.50e-04 8.28e-04 8.3e-02 1.4e+00 9.9e-02 1.84e-03
## 11 22 -4.151e+03 2.83e-03 2.97e-03 6.5e-02 0.0e+00 9.7e-02 2.97e-03
## 12 24 -4.153e+03 4.75e-04 1.83e-03 3.3e-02 1.3e+00 4.8e-02 2.46e-03
## 13 25 -4.154e+03 1.57e-04 1.14e-04 2.9e-03 0.0e+00 5.0e-03 1.14e-04
## 14 28 -4.154e+03 7.58e-05 9.14e-05 1.8e-02 2.9e-01 2.6e-02 9.53e-05
## 15 29 -4.154e+03 2.44e-05 1.44e-05 3.9e-03 0.0e+00 7.5e-03 1.44e-05
## 16 30 -4.154e+03 3.02e-05 3.32e-05 1.2e-02 0.0e+00 2.3e-02 3.32e-05
## 17 31 -4.154e+03 5.38e-07 3.86e-07 1.3e-03 0.0e+00 2.7e-03 3.86e-07
## 18 50 -4.154e+03 -4.60e-15 5.44e-15 1.1e-14 1.5e+06 1.6e-14 1.08e-08
##
## ***** FALSE CONVERGENCE *****
##
## FUNCTION -4.154077e+03 RELDX 1.082e-14
## FUNC. EVALS 50 GRAD. EVALS 18
## PRELDF 5.437e-15 NPRELDF 1.079e-08
##
## I FINAL X(I) D(I) G(I)
##
## 1 4.899892e-05 1.000e+00 1.421e+03
## 2 1.762710e-01 1.000e+00 4.281e-01
## 3 7.350063e-01 1.000e+00 4.144e-01fit$coe[3]## b1
## 0.7350063#so take the value of lambda as 0.74
#S&P500 index price
lp.GSPC = log(GSPC[,5])
r11 = diff(lp.GSPC)
r11 = r11[2:length(r11)]
fit = garch(r11)##
## ***** ESTIMATION WITH ANALYTICAL GRADIENT *****
##
##
## I INITIAL X(I) D(I)
##
## 1 1.733221e-04 1.000e+00
## 2 5.000000e-02 1.000e+00
## 3 5.000000e-02 1.000e+00
##
## IT NF F RELDF PRELDF RELDX STPPAR D*STEP NPRELDF
## 0 1 -4.751e+03
## 1 7 -4.762e+03 2.22e-03 2.73e-03 2.0e-04 3.3e+10 2.0e-05 4.53e+07
## 2 8 -4.765e+03 7.59e-04 1.06e-03 2.0e-04 2.4e+00 2.0e-05 1.11e+02
## 3 9 -4.765e+03 1.77e-05 2.25e-05 1.9e-04 2.0e+00 2.0e-05 1.02e+02
## 4 16 -4.810e+03 9.19e-03 1.39e-02 4.0e-01 2.0e+00 6.8e-02 1.02e+02
## 5 17 -4.824e+03 2.96e-03 4.30e-03 2.8e-01 2.0e+00 6.8e-02 6.81e+00
## 6 20 -4.885e+03 1.25e-02 1.41e-02 5.1e-01 2.0e+00 2.6e-01 7.43e+00
## 7 21 -4.966e+03 1.64e-02 1.73e-02 2.5e-01 2.0e+00 2.6e-01 1.09e+00
## 8 23 -4.982e+03 3.17e-03 4.02e-03 2.0e-02 2.0e+00 2.6e-02 4.53e+00
## 9 24 -4.992e+03 1.92e-03 2.64e-03 1.9e-02 2.0e+00 2.6e-02 2.93e+00
## 10 26 -5.004e+03 2.40e-03 4.72e-03 5.5e-02 2.0e+00 8.1e-02 4.10e+00
## 11 36 -5.005e+03 3.12e-04 2.26e-03 1.8e-06 3.1e+00 2.8e-06 6.53e-02
## 12 37 -5.007e+03 3.52e-04 2.58e-04 1.7e-06 2.0e+00 2.8e-06 2.84e-03
## 13 38 -5.007e+03 3.09e-05 4.31e-05 1.8e-06 2.0e+00 2.8e-06 2.93e-04
## 14 39 -5.007e+03 2.96e-06 2.96e-06 1.8e-06 2.0e+00 2.8e-06 4.94e-06
## 15 49 -5.007e+03 2.86e-09 1.22e-08 3.0e-09 2.6e+00 4.6e-09 5.74e-06
## 16 62 -5.007e+03 -1.62e-14 3.91e-14 1.7e-14 1.1e+05 2.6e-14 5.82e-06
##
## ***** FALSE CONVERGENCE *****
##
## FUNCTION -5.007264e+03 RELDX 1.699e-14
## FUNC. EVALS 62 GRAD. EVALS 16
## PRELDF 3.907e-14 NPRELDF 5.823e-06
##
## I FINAL X(I) D(I) G(I)
##
## 1 5.299235e-06 1.000e+00 7.507e+03
## 2 2.144131e-01 1.000e+00 1.283e+00
## 3 7.669003e-01 1.000e+00 -1.199e+00fit$coe[3]## b1
## 0.7669003#so take the value of lambda as 0.79
#λ for period 2013-2017
#Apple Inc
lp.AAPL2 = log(AAPL2[,5])
r1 = diff(lp.AAPL2)
r1 = r1[2:length(r1)]
fit = garch(r1)##
## ***** ESTIMATION WITH ANALYTICAL GRADIENT *****
##
##
## I INITIAL X(I) D(I)
##
## 1 2.070021e-04 1.000e+00
## 2 5.000000e-02 1.000e+00
## 3 5.000000e-02 1.000e+00
##
## IT NF F RELDF PRELDF RELDX STPPAR D*STEP NPRELDF
## 0 1 -4.568e+03
## 1 8 -4.568e+03 2.05e-05 3.95e-05 3.7e-05 1.3e+10 3.7e-06 2.66e+05
## 2 17 -4.574e+03 1.27e-03 2.24e-03 4.8e-01 2.0e+00 9.2e-02 1.35e+00
## 3 20 -4.580e+03 1.23e-03 1.34e-03 7.7e-01 1.6e+00 3.7e-01 2.95e-02
## 4 22 -4.582e+03 4.47e-04 3.84e-04 8.0e-02 2.0e+00 7.4e-02 1.46e-01
## 5 23 -4.584e+03 5.63e-04 7.80e-04 1.3e-01 2.0e+00 1.5e-01 1.28e+00
## 6 25 -4.585e+03 9.70e-05 3.25e-04 1.1e-02 1.9e+00 1.5e-02 9.33e-04
## 7 26 -4.586e+03 1.97e-04 2.03e-04 1.0e-02 1.7e+00 1.5e-02 4.67e-04
## 8 29 -4.589e+03 7.04e-04 4.68e-04 6.9e-02 0.0e+00 1.2e-01 4.68e-04
## 9 31 -4.589e+03 4.92e-05 2.72e-04 1.5e-02 1.7e+00 2.4e-02 1.96e-03
## 10 32 -4.590e+03 5.45e-05 7.11e-05 1.5e-02 7.9e-01 2.4e-02 1.04e-04
## 11 34 -4.590e+03 3.41e-07 6.97e-05 4.8e-03 1.5e+00 7.7e-03 9.12e-05
## 12 35 -4.590e+03 2.99e-05 3.04e-05 1.9e-03 1.5e+00 3.9e-03 3.10e-05
## 13 36 -4.590e+03 1.44e-07 1.12e-07 6.8e-04 0.0e+00 1.5e-03 1.12e-07
## 14 58 -4.590e+03 -5.35e-15 8.69e-15 7.3e-15 8.2e+05 1.2e-14 3.72e-08
##
## ***** FALSE CONVERGENCE *****
##
## FUNCTION -4.589667e+03 RELDX 7.275e-15
## FUNC. EVALS 58 GRAD. EVALS 14
## PRELDF 8.685e-15 NPRELDF 3.724e-08
##
## I FINAL X(I) D(I) G(I)
##
## 1 2.563826e-05 1.000e+00 -3.440e+03
## 2 9.732144e-02 1.000e+00 -1.325e+00
## 3 7.964072e-01 1.000e+00 -8.278e-01fit$coe[3]## b1
## 0.7964072#so take the value of lambda as 0.80
#Amazon.com Inc
lp.AMZN2 = log(AMZN2[,5])
r2 = diff(lp.AMZN2)
r2 = r2[2:length(r2)]
fit = garch(r2)##
## ***** ESTIMATION WITH ANALYTICAL GRADIENT *****
##
##
## I INITIAL X(I) D(I)
##
## 1 3.023619e-04 1.000e+00
## 2 5.000000e-02 1.000e+00
## 3 5.000000e-02 1.000e+00
##
## IT NF F RELDF PRELDF RELDX STPPAR D*STEP NPRELDF
## 0 1 -4.338e+03
## 1 7 -4.338e+03 5.19e-05 1.21e-04 1.0e-04 5.2e+09 1.0e-05 3.18e+05
## 2 8 -4.338e+03 1.86e-06 2.04e-06 9.9e-05 2.0e+00 1.0e-05 4.02e+00
## 3 15 -4.344e+03 1.32e-03 1.76e-03 2.9e-01 2.0e+00 4.1e-02 4.02e+00
## 4 19 -4.349e+03 1.12e-03 5.44e-04 8.4e-01 4.4e-01 5.3e-01 1.83e-03
## 5 27 -4.349e+03 2.16e-05 1.83e-04 5.6e-06 4.4e+00 6.5e-06 1.62e+00
## 6 28 -4.349e+03 3.52e-05 3.16e-05 5.1e-06 2.0e+00 6.5e-06 1.78e+00
## 7 37 -4.350e+03 3.68e-04 3.46e-04 6.0e-02 2.0e+00 7.4e-02 1.97e+00
## 8 39 -4.351e+03 1.89e-04 2.08e-04 3.3e-02 2.0e+00 4.4e-02 2.49e+00
## 9 45 -4.351e+03 2.74e-07 3.77e-05 1.6e-06 5.1e+00 2.2e-06 2.33e-03
## 10 46 -4.351e+03 1.01e-05 9.33e-06 7.5e-07 2.0e+00 1.1e-06 9.28e-04
## 11 55 -4.352e+03 5.81e-05 6.50e-05 1.3e-02 1.3e+00 1.8e-02 9.57e-04
## 12 56 -4.352e+03 6.18e-05 5.82e-05 1.2e-02 5.0e-01 1.8e-02 6.91e-05
## 13 57 -4.352e+03 6.29e-05 4.09e-05 6.0e-03 0.0e+00 8.8e-03 4.09e-05
## 14 58 -4.352e+03 5.16e-05 4.70e-05 1.0e-02 2.6e-01 1.9e-02 4.92e-05
## 15 59 -4.352e+03 8.71e-06 8.61e-06 7.0e-03 0.0e+00 1.2e-02 8.61e-06
## 16 60 -4.352e+03 6.91e-09 1.50e-07 1.1e-03 0.0e+00 1.7e-03 1.50e-07
## 17 63 -4.352e+03 9.12e-10 9.12e-11 1.9e-06 1.9e+00 2.9e-06 8.29e-10
## 18 72 -4.352e+03 -1.50e-14 8.05e-17 2.1e-14 4.9e+03 2.9e-14 7.41e-10
##
## ***** FALSE CONVERGENCE *****
##
## FUNCTION -4.352387e+03 RELDX 2.062e-14
## FUNC. EVALS 72 GRAD. EVALS 18
## PRELDF 8.049e-17 NPRELDF 7.412e-10
##
## I FINAL X(I) D(I) G(I)
##
## 1 5.874730e-05 1.000e+00 -1.199e+01
## 2 1.359097e-01 1.000e+00 -1.341e-02
## 3 7.086053e-01 1.000e+00 4.308e-02fit$coe[3]## b1
## 0.7086053#so take the value of lambda as 0.71
#Microsoft Corp
lp.MSFT2 = log(MSFT2[,5])
r3 = diff(lp.MSFT2)
r3 = r3[2:length(r3)]
fit = garch(r3)##
## ***** ESTIMATION WITH ANALYTICAL GRADIENT *****
##
##
## I INITIAL X(I) D(I)
##
## 1 1.791948e-04 1.000e+00
## 2 5.000000e-02 1.000e+00
## 3 5.000000e-02 1.000e+00
##
## IT NF F RELDF PRELDF RELDX STPPAR D*STEP NPRELDF
## 0 1 -4.655e+03
## 1 8 -4.655e+03 7.13e-06 1.35e-05 1.8e-05 1.9e+10 1.8e-06 1.25e+05
## 2 9 -4.655e+03 5.44e-08 5.77e-08 1.8e-05 2.0e+00 1.8e-06 9.43e-01
## 3 17 -4.659e+03 8.26e-04 1.21e-03 3.8e-01 2.0e+00 6.0e-02 9.42e-01
## 4 20 -4.663e+03 8.82e-04 9.75e-04 6.8e-01 1.8e+00 2.4e-01 3.94e-02
## 5 22 -4.666e+03 5.32e-04 5.70e-04 2.0e-01 1.7e+00 1.5e-01 4.96e-03
## 6 24 -4.668e+03 5.20e-04 5.82e-04 1.5e-01 1.6e+00 1.5e-01 1.07e-02
## 7 30 -4.668e+03 3.17e-05 5.36e-05 1.2e-06 1.4e+01 1.4e-06 4.77e-04
## 8 31 -4.668e+03 6.18e-07 8.24e-07 1.1e-06 2.0e+00 1.4e-06 4.65e-05
## 9 39 -4.668e+03 9.30e-06 1.40e-05 9.6e-03 1.1e+00 1.2e-02 4.37e-05
## 10 40 -4.668e+03 6.40e-06 5.25e-06 5.0e-03 0.0e+00 6.6e-03 5.25e-06
## 11 42 -4.668e+03 4.07e-06 3.01e-06 6.0e-03 0.0e+00 9.5e-03 3.01e-06
## 12 43 -4.668e+03 8.56e-07 5.01e-07 3.7e-03 0.0e+00 5.4e-03 5.01e-07
## 13 44 -4.668e+03 3.90e-08 2.81e-09 2.8e-04 0.0e+00 3.5e-04 2.81e-09
## 14 45 -4.668e+03 2.82e-09 3.42e-11 1.4e-05 0.0e+00 2.2e-05 3.42e-11
## 15 46 -4.668e+03 2.43e-11 2.18e-13 1.8e-06 0.0e+00 2.4e-06 2.18e-13
## 16 47 -4.668e+03 -8.29e-12 2.35e-16 8.0e-08 0.0e+00 1.1e-07 2.35e-16
##
## ***** RELATIVE FUNCTION CONVERGENCE *****
##
## FUNCTION -4.668292e+03 RELDX 8.039e-08
## FUNC. EVALS 47 GRAD. EVALS 16
## PRELDF 2.347e-16 NPRELDF 2.347e-16
##
## I FINAL X(I) D(I) G(I)
##
## 1 4.971301e-05 1.000e+00 1.525e-01
## 2 1.798582e-01 1.000e+00 2.633e-05
## 3 6.043928e-01 1.000e+00 1.364e-05fit$coe[3]## b1
## 0.6043928#so take the value of lambda as 0.60
#Alphabet Inc
lp.GOOG2 = log(GOOG2[,5])
r4 = diff(lp.GOOG2)
r4 = r4[2:length(r4)]
fit = garch(r4)##
## ***** ESTIMATION WITH ANALYTICAL GRADIENT *****
##
##
## I INITIAL X(I) D(I)
##
## 1 1.714040e-04 1.000e+00
## 2 5.000000e-02 1.000e+00
## 3 5.000000e-02 1.000e+00
##
## IT NF F RELDF PRELDF RELDX STPPAR D*STEP NPRELDF
## 0 1 -4.692e+03
## 1 7 -4.693e+03 9.63e-05 2.95e-04 1.0e-04 1.4e+10 1.0e-05 2.04e+06
## 2 8 -4.693e+03 1.57e-05 1.73e-05 9.6e-05 2.0e+00 1.0e-05 8.76e+00
## 3 16 -4.713e+03 4.18e-03 7.37e-03 5.5e-01 2.0e+00 1.2e-01 8.74e+00
## 4 19 -4.736e+03 4.84e-03 3.77e-03 8.2e-01 1.7e+00 4.8e-01 1.20e-01
## 5 28 -4.738e+03 5.53e-04 1.50e-03 8.1e-06 4.6e+00 8.7e-06 8.63e-02
## 6 29 -4.739e+03 5.52e-05 4.61e-05 7.9e-06 2.0e+00 8.7e-06 2.46e-02
## 7 30 -4.739e+03 2.03e-06 2.40e-06 8.1e-06 2.0e+00 8.7e-06 3.28e-02
## 8 38 -4.740e+03 2.80e-04 6.49e-04 8.2e-02 1.9e+00 9.6e-02 3.23e-02
## 9 39 -4.740e+03 3.96e-05 4.09e-05 2.5e-02 0.0e+00 3.1e-02 4.09e-05
## 10 40 -4.740e+03 1.41e-05 9.30e-06 5.8e-03 0.0e+00 8.9e-03 9.30e-06
## 11 41 -4.740e+03 6.85e-06 6.95e-06 8.1e-03 0.0e+00 1.0e-02 6.95e-06
## 12 42 -4.740e+03 1.43e-07 6.70e-08 1.2e-03 0.0e+00 1.5e-03 6.70e-08
## 13 43 -4.740e+03 1.94e-09 1.43e-09 4.2e-05 0.0e+00 4.9e-05 1.43e-09
## 14 44 -4.740e+03 8.39e-10 9.79e-12 6.8e-06 0.0e+00 8.9e-06 9.79e-12
## 15 45 -4.740e+03 3.99e-11 9.18e-15 4.6e-07 0.0e+00 5.8e-07 9.18e-15
## 16 46 -4.740e+03 7.11e-13 4.49e-18 9.6e-09 0.0e+00 1.2e-08 4.49e-18
## 17 47 -4.740e+03 -5.76e-16 3.21e-22 2.4e-11 0.0e+00 3.3e-11 3.21e-22
##
## ***** X- AND RELATIVE FUNCTION CONVERGENCE *****
##
## FUNCTION -4.740151e+03 RELDX 2.404e-11
## FUNC. EVALS 47 GRAD. EVALS 17
## PRELDF 3.210e-22 NPRELDF 3.210e-22
##
## I FINAL X(I) D(I) G(I)
##
## 1 3.801146e-05 1.000e+00 5.330e-04
## 2 2.464730e-01 1.000e+00 7.904e-08
## 3 5.896038e-01 1.000e+00 1.040e-07fit$coe[3]## b1
## 0.5896038#so take the value of lambda as 0.59
#Tesla Inc
lp.TSLA2 = log(TSLA2[,5])
r5 = diff(lp.TSLA2)
r5 = r5[2:length(r5)]
fit = garch(r5)##
## ***** ESTIMATION WITH ANALYTICAL GRADIENT *****
##
##
## I INITIAL X(I) D(I)
##
## 1 7.935678e-04 1.000e+00
## 2 5.000000e-02 1.000e+00
## 3 5.000000e-02 1.000e+00
##
## IT NF F RELDF PRELDF RELDX STPPAR D*STEP NPRELDF
## 0 1 -3.741e+03
## 1 7 -3.741e+03 4.29e-05 7.56e-05 1.7e-04 1.0e+09 1.7e-05 3.84e+04
## 2 8 -3.741e+03 1.15e-06 1.20e-06 1.6e-04 2.0e+00 1.7e-05 1.93e+00
## 3 15 -3.745e+03 9.67e-04 1.60e-03 3.3e-01 2.0e+00 5.0e-02 1.92e+00
## 4 18 -3.756e+03 2.92e-03 2.57e-03 6.2e-01 1.9e+00 2.0e-01 3.72e-01
## 5 20 -3.764e+03 2.19e-03 2.60e-03 2.4e-01 2.0e+00 1.7e-01 2.69e+01
## 6 29 -3.764e+03 7.01e-05 1.38e-04 1.7e-05 5.4e+00 1.4e-05 4.00e-03
## 7 30 -3.764e+03 6.84e-07 7.84e-07 1.7e-05 2.0e+00 1.4e-05 2.49e-03
## 8 38 -3.767e+03 8.37e-04 1.48e-03 1.7e-01 8.0e-01 1.8e-01 2.50e-03
## 9 39 -3.768e+03 2.69e-04 2.15e-04 2.0e-02 0.0e+00 3.0e-02 2.15e-04
## 10 42 -3.772e+03 9.88e-04 4.41e-04 8.4e-02 0.0e+00 1.3e-01 4.41e-04
## 11 44 -3.775e+03 7.18e-04 6.50e-04 3.1e-02 2.0e+00 5.1e-02 4.24e-02
## 12 46 -3.777e+03 6.11e-04 6.16e-04 2.9e-02 2.0e+00 5.1e-02 1.88e-01
## 13 48 -3.778e+03 2.56e-04 2.23e-04 8.7e-03 2.0e+00 1.6e-02 7.46e-01
## 14 50 -3.781e+03 7.45e-04 5.73e-04 1.7e-02 2.0e+00 3.2e-02 5.33e+00
## 15 52 -3.789e+03 2.21e-03 1.69e-03 3.2e-02 2.0e+00 6.3e-02 4.13e+02
## 16 54 -3.791e+03 5.51e-04 5.29e-04 6.2e-03 2.0e+00 1.3e-02 5.23e+04
## 17 55 -3.792e+03 2.81e-04 1.91e-03 1.2e-02 2.0e+00 2.5e-02 3.56e+04
## 18 58 -3.793e+03 9.52e-05 1.56e-03 5.9e-04 2.0e+00 1.2e-03 1.24e+03
## 19 59 -3.796e+03 7.97e-04 1.29e-03 2.9e-04 2.0e+00 6.1e-04 3.87e+03
## 20 60 -3.796e+03 6.16e-05 5.56e-05 2.9e-04 2.0e+00 6.1e-04 8.67e+00
## 21 61 -3.796e+03 1.07e-05 1.27e-05 2.9e-04 2.0e+00 6.1e-04 5.30e+02
## 22 63 -3.796e+03 1.27e-05 2.03e-05 7.4e-04 2.0e+00 1.6e-03 1.74e+02
## 23 71 -3.796e+03 1.09e-06 1.89e-06 4.3e-08 4.9e+00 8.4e-08 1.03e+00
## 24 72 -3.796e+03 4.39e-10 5.72e-10 4.0e-08 2.0e+00 8.4e-08 2.35e-01
## 25 79 -3.796e+03 2.08e-06 3.39e-06 3.3e-04 2.0e+00 6.9e-04 2.35e-01
## 26 81 -3.796e+03 1.30e-05 2.58e-05 1.1e-03 1.8e+00 2.7e-03 3.87e-04
## 27 82 -3.796e+03 9.28e-08 1.10e-07 2.4e-05 0.0e+00 6.4e-05 1.10e-07
## 28 83 -3.796e+03 4.13e-09 5.10e-09 2.1e-06 0.0e+00 5.0e-06 5.10e-09
## 29 84 -3.796e+03 1.36e-10 4.05e-10 1.1e-06 0.0e+00 2.4e-06 4.05e-10
## 30 85 -3.796e+03 5.56e-11 1.25e-11 3.6e-07 0.0e+00 7.6e-07 1.25e-11
## 31 86 -3.796e+03 -1.66e-12 5.17e-14 2.7e-08 0.0e+00 5.2e-08 5.17e-14
##
## ***** RELATIVE FUNCTION CONVERGENCE *****
##
## FUNCTION -3.796244e+03 RELDX 2.695e-08
## FUNC. EVALS 86 GRAD. EVALS 31
## PRELDF 5.165e-14 NPRELDF 5.165e-14
##
## I FINAL X(I) D(I) G(I)
##
## 1 5.011792e-06 1.000e+00 1.247e+01
## 2 2.514991e-02 1.000e+00 1.240e-02
## 3 9.690119e-01 1.000e+00 1.313e-02fit$coe[3]## b1
## 0.9690119#so take the value of lambda as 0.97
#NVIDIA Corp
lp.NVDA2 = log(NVDA2[,5])
r6 = diff(lp.NVDA2)
r6 = r6[2:length(r6)]
fit = garch(r6)##
## ***** ESTIMATION WITH ANALYTICAL GRADIENT *****
##
##
## I INITIAL X(I) D(I)
##
## 1 4.164192e-04 1.000e+00
## 2 5.000000e-02 1.000e+00
## 3 5.000000e-02 1.000e+00
##
## IT NF F RELDF PRELDF RELDX STPPAR D*STEP NPRELDF
## 0 1 -4.133e+03
## 1 8 -4.133e+03 6.05e-06 1.18e-05 3.9e-05 3.3e+09 3.9e-06 1.94e+04
## 2 16 -4.139e+03 1.39e-03 1.91e-03 3.9e-01 2.0e+00 6.3e-02 1.90e+00
## 3 19 -4.142e+03 7.36e-04 7.69e-04 7.1e-01 1.6e+00 2.5e-01 1.87e-02
## 4 21 -4.142e+03 1.23e-04 2.53e-04 9.3e-02 1.3e+00 6.4e-02 5.39e-04
## 5 22 -4.143e+03 1.09e-04 2.15e-04 5.6e-02 0.0e+00 5.6e-02 2.15e-04
## 6 24 -4.144e+03 2.59e-04 2.48e-04 2.5e-01 4.4e-01 1.7e-01 2.70e-04
## 7 25 -4.144e+03 4.67e-05 9.01e-05 2.8e-01 0.0e+00 1.3e-01 9.01e-05
## 8 26 -4.144e+03 6.99e-05 6.62e-05 1.3e-01 0.0e+00 6.2e-02 6.62e-05
## 9 27 -4.144e+03 1.03e-05 7.29e-06 2.4e-02 0.0e+00 1.3e-02 7.29e-06
## 10 28 -4.144e+03 3.06e-06 2.84e-06 1.9e-02 0.0e+00 1.3e-02 2.84e-06
## 11 29 -4.144e+03 9.03e-08 7.17e-08 1.4e-03 0.0e+00 1.0e-03 7.17e-08
## 12 30 -4.144e+03 8.37e-09 6.15e-09 4.5e-04 0.0e+00 2.4e-04 6.15e-09
## 13 31 -4.144e+03 3.35e-10 1.64e-10 1.1e-04 0.0e+00 7.1e-05 1.64e-10
## 14 32 -4.144e+03 -3.54e-12 1.49e-12 1.3e-05 0.0e+00 6.9e-06 1.49e-12
##
## ***** RELATIVE FUNCTION CONVERGENCE *****
##
## FUNCTION -4.144133e+03 RELDX 1.288e-05
## FUNC. EVALS 32 GRAD. EVALS 14
## PRELDF 1.493e-12 NPRELDF 1.493e-12
##
## I FINAL X(I) D(I) G(I)
##
## 1 2.953302e-04 1.000e+00 5.879e-01
## 2 2.643383e-01 1.000e+00 2.015e-03
## 3 1.705457e-01 1.000e+00 8.064e-04fit$coe[3]## b1
## 0.1705457#so take the value of lambda as 0.17
#PepsiCo Inc
lp.PEP2 = log(PEP2[,5])
r7 = diff(lp.PEP2)
r7 = r7[2:length(r7)]
fit = garch(r7)##
## ***** ESTIMATION WITH ANALYTICAL GRADIENT *****
##
##
## I INITIAL X(I) D(I)
##
## 1 6.242552e-05 1.000e+00
## 2 5.000000e-02 1.000e+00
## 3 5.000000e-02 1.000e+00
##
## IT NF F RELDF PRELDF RELDX STPPAR D*STEP NPRELDF
## 0 1 -5.314e+03
## 1 8 -5.314e+03 7.05e-05 1.21e-04 2.0e-05 1.6e+11 2.0e-06 9.94e+06
## 2 9 -5.314e+03 1.52e-06 1.66e-06 1.9e-05 2.0e+00 2.0e-06 4.08e+00
## 3 18 -5.324e+03 1.94e-03 3.37e-03 4.6e-01 2.0e+00 8.6e-02 4.07e+00
## 4 21 -5.330e+03 1.01e-03 1.22e-03 6.9e-01 1.9e+00 2.5e-01 5.80e-02
## 5 28 -5.330e+03 2.13e-05 4.31e-05 1.4e-06 2.3e+01 8.6e-07 7.71e-04
## 6 29 -5.330e+03 2.19e-08 4.34e-08 1.4e-06 2.0e+00 8.6e-07 6.08e-04
## 7 30 -5.330e+03 1.89e-09 2.87e-09 1.4e-06 2.0e+00 8.6e-07 6.09e-04
## 8 39 -5.331e+03 2.69e-04 3.12e-04 1.6e-01 8.7e-01 1.1e-01 6.09e-04
## 9 40 -5.331e+03 1.86e-05 6.09e-05 8.0e-02 0.0e+00 7.3e-02 6.09e-05
## 10 41 -5.332e+03 2.69e-05 2.57e-05 4.7e-02 0.0e+00 4.4e-02 2.57e-05
## 11 42 -5.332e+03 2.90e-06 2.06e-06 5.7e-03 0.0e+00 6.0e-03 2.06e-06
## 12 43 -5.332e+03 1.41e-06 1.30e-06 5.2e-03 1.3e-01 6.0e-03 1.32e-06
## 13 44 -5.332e+03 5.23e-08 1.13e-08 1.1e-03 0.0e+00 9.9e-04 1.13e-08
## 14 45 -5.332e+03 3.36e-09 4.77e-11 4.8e-05 0.0e+00 4.2e-05 4.77e-11
## 15 46 -5.332e+03 4.05e-10 6.22e-13 5.6e-06 0.0e+00 5.0e-06 6.22e-13
## 16 47 -5.332e+03 2.90e-12 2.79e-16 5.0e-08 0.0e+00 5.3e-08 2.79e-16
## 17 48 -5.332e+03 -4.47e-14 2.64e-19 1.5e-09 0.0e+00 1.5e-09 2.64e-19
##
## ***** X- AND RELATIVE FUNCTION CONVERGENCE *****
##
## FUNCTION -5.331560e+03 RELDX 1.502e-09
## FUNC. EVALS 48 GRAD. EVALS 17
## PRELDF 2.640e-19 NPRELDF 2.640e-19
##
## I FINAL X(I) D(I) G(I)
##
## 1 2.683234e-05 1.000e+00 2.654e-02
## 2 1.799625e-01 1.000e+00 2.243e-06
## 3 4.345448e-01 1.000e+00 1.865e-06fit$coe[3]## b1
## 0.4345448#so take the value of lambda as 0.43
#Meta Platforms Inc
lp.META2 = log(META2[,5])
r8 = diff(lp.META2)
r8 = r8[2:length(r8)]
fit = garch(r8)##
## ***** ESTIMATION WITH ANALYTICAL GRADIENT *****
##
##
## I INITIAL X(I) D(I)
##
## 1 3.547507e-04 1.000e+00
## 2 5.000000e-02 1.000e+00
## 3 5.000000e-02 1.000e+00
##
## IT NF F RELDF PRELDF RELDX STPPAR D*STEP NPRELDF
## 0 1 -4.232e+03
## 1 8 -4.232e+03 6.82e-06 1.33e-05 3.5e-05 4.5e+09 3.5e-06 3.02e+04
## 2 17 -4.242e+03 2.22e-03 3.94e-03 5.8e-01 2.0e+00 1.4e-01 1.72e+00
## 3 18 -4.243e+03 4.50e-04 5.23e-04 3.7e-01 1.9e+00 1.4e-01 1.53e-02
## 4 20 -4.244e+03 1.79e-04 2.37e-04 2.5e-01 1.8e+00 1.3e-01 3.16e-03
## 5 22 -4.246e+03 4.65e-04 2.87e-04 1.7e-01 8.9e-01 1.3e-01 1.63e-03
## 6 23 -4.247e+03 2.27e-04 1.18e-03 2.2e-01 2.0e+00 2.7e-01 4.16e-01
## 7 32 -4.251e+03 8.09e-04 1.78e-03 6.4e-06 5.0e+00 9.5e-06 7.04e-03
## 8 33 -4.251e+03 1.22e-05 9.82e-06 6.3e-06 2.0e+00 9.5e-06 2.85e-03
## 9 34 -4.251e+03 7.21e-07 6.77e-07 6.3e-06 2.0e+00 9.5e-06 2.19e-03
## 10 41 -4.251e+03 1.18e-04 2.03e-04 2.1e-02 1.8e+00 3.0e-02 2.21e-03
## 11 44 -4.255e+03 9.82e-04 4.36e-04 6.6e-02 0.0e+00 1.3e-01 4.36e-04
## 12 46 -4.258e+03 5.42e-04 4.57e-04 2.5e-02 2.0e+00 5.1e-02 3.56e-02
## 13 47 -4.261e+03 7.77e-04 1.20e-03 4.7e-02 2.0e+00 1.0e-01 4.23e+00
## 14 55 -4.262e+03 3.61e-04 7.71e-04 8.6e-07 8.0e+00 1.6e-06 5.90e-03
## 15 56 -4.262e+03 2.19e-06 1.75e-06 8.3e-07 2.0e+00 1.6e-06 1.02e-02
## 16 65 -4.273e+03 2.53e-03 1.58e-03 1.3e-02 8.6e-01 2.6e-02 1.05e-02
## 17 67 -4.276e+03 6.09e-04 6.05e-04 2.6e-03 2.0e+00 5.3e-03 1.66e+00
## 18 69 -4.282e+03 1.35e-03 1.24e-03 5.1e-03 2.0e+00 1.1e-02 1.24e+01
## 19 75 -4.283e+03 2.51e-04 3.48e-04 1.8e-07 9.2e+00 3.4e-07 2.95e+01
## 20 76 -4.283e+03 2.57e-05 2.88e-05 1.4e-07 2.0e+00 3.4e-07 6.43e+01
## 21 77 -4.283e+03 5.50e-07 7.27e-07 1.6e-07 2.0e+00 3.4e-07 5.86e+01
## 22 86 -4.287e+03 1.04e-03 9.86e-04 2.7e-03 2.0e+00 5.6e-03 5.88e+01
## 23 88 -4.288e+03 2.09e-04 2.08e-04 5.4e-04 2.0e+00 1.1e-03 2.31e+03
## 24 94 -4.288e+03 4.50e-07 2.01e-06 1.1e-08 5.6e+01 2.2e-08 3.86e+00
## 25 105 -4.292e+03 9.80e-04 1.19e-03 2.8e-03 2.0e+00 5.9e-03 3.29e+00
## 26 106 -4.298e+03 1.33e-03 2.47e-03 2.5e-03 1.9e+00 5.9e-03 8.05e-02
## 27 113 -4.298e+03 2.64e-05 6.97e-05 4.8e-08 8.4e+00 9.5e-08 4.30e-04
## 28 114 -4.298e+03 1.99e-06 1.70e-06 3.7e-08 2.0e+00 9.5e-08 1.37e-04
## 29 122 -4.299e+03 6.24e-05 7.47e-05 6.2e-04 8.2e-01 1.6e-03 1.31e-04
## 30 123 -4.299e+03 1.54e-04 1.39e-04 1.3e-03 0.0e+00 3.1e-03 1.39e-04
## 31 124 -4.299e+03 7.64e-06 5.94e-06 3.1e-04 0.0e+00 7.7e-04 5.94e-06
## 32 125 -4.299e+03 9.24e-07 8.42e-07 1.3e-04 0.0e+00 3.0e-04 8.42e-07
## 33 126 -4.299e+03 1.11e-08 1.24e-08 1.1e-05 0.0e+00 2.4e-05 1.24e-08
## 34 135 -4.299e+03 -2.12e-13 6.59e-14 2.1e-14 5.0e+03 4.2e-14 1.75e-10
##
## ***** FALSE CONVERGENCE *****
##
## FUNCTION -4.299216e+03 RELDX 2.123e-14
## FUNC. EVALS 135 GRAD. EVALS 34
## PRELDF 6.592e-14 NPRELDF 1.748e-10
##
## I FINAL X(I) D(I) G(I)
##
## 1 5.705821e-07 1.000e+00 -6.794e+03
## 2 1.661282e-02 1.000e+00 -1.395e+00
## 3 9.823599e-01 1.000e+00 -1.908e+00fit$coe[3]## b1
## 0.9823599#so take the value of lambda as 0.98
#Costco Wholesale Corp
lp.COST2 = log(COST2[,5])
r9 = diff(lp.COST2)
r9 = r9[2:length(r9)]
fit = garch(r9)##
## ***** ESTIMATION WITH ANALYTICAL GRADIENT *****
##
##
## I INITIAL X(I) D(I)
##
## 1 1.016876e-04 1.000e+00
## 2 5.000000e-02 1.000e+00
## 3 5.000000e-02 1.000e+00
##
## IT NF F RELDF PRELDF RELDX STPPAR D*STEP NPRELDF
## 0 1 -5.005e+03
## 1 8 -5.005e+03 7.42e-06 1.34e-05 1.1e-05 6.1e+10 1.1e-06 4.05e+05
## 2 9 -5.005e+03 9.99e-08 1.05e-07 1.0e-05 2.0e+00 1.1e-06 8.76e-01
## 3 18 -5.009e+03 8.72e-04 1.29e-03 4.1e-01 2.0e+00 6.9e-02 8.75e-01
## 4 20 -5.013e+03 7.15e-04 1.73e-03 6.9e-01 2.0e+00 2.8e-01 1.11e-01
## 5 27 -5.013e+03 5.82e-05 3.42e-04 8.2e-06 3.2e+00 5.4e-06 4.91e-04
## 6 28 -5.013e+03 6.33e-05 5.74e-05 7.5e-06 2.0e+00 5.4e-06 1.07e-04
## 7 29 -5.013e+03 4.73e-07 5.76e-07 8.1e-06 2.0e+00 5.4e-06 3.25e-05
## 8 36 -5.013e+03 1.88e-05 2.09e-05 3.4e-02 7.6e-01 2.2e-02 3.30e-05
## 9 37 -5.013e+03 9.05e-06 6.38e-06 3.6e-02 0.0e+00 2.2e-02 6.38e-06
## 10 38 -5.013e+03 3.39e-06 2.46e-06 2.3e-02 0.0e+00 1.4e-02 2.46e-06
## 11 39 -5.013e+03 8.53e-07 6.65e-07 6.7e-03 0.0e+00 4.7e-03 6.65e-07
## 12 40 -5.013e+03 8.34e-08 1.06e-07 2.1e-03 0.0e+00 1.4e-03 1.06e-07
## 13 53 -5.013e+03 -1.45e-15 1.01e-16 7.0e-15 2.9e+05 3.8e-15 2.27e-09
##
## ***** FALSE CONVERGENCE *****
##
## FUNCTION -5.013443e+03 RELDX 6.994e-15
## FUNC. EVALS 53 GRAD. EVALS 13
## PRELDF 1.011e-16 NPRELDF 2.268e-09
##
## I FINAL X(I) D(I) G(I)
##
## 1 6.428952e-05 1.000e+00 1.336e+02
## 2 1.847318e-01 1.000e+00 -4.745e-02
## 3 2.713638e-01 1.000e+00 -3.851e-02fit$coe[3]## b1
## 0.2713638#so take the value of lambda as 0.27
#Broadcom Inc
lp.AVGO2 = log(AVGO2[,5])
r10 = diff(lp.AVGO2)
r10 = r10[2:length(r10)]
fit = garch(r10)##
## ***** ESTIMATION WITH ANALYTICAL GRADIENT *****
##
##
## I INITIAL X(I) D(I)
##
## 1 3.570765e-04 1.000e+00
## 2 5.000000e-02 1.000e+00
## 3 5.000000e-02 1.000e+00
##
## IT NF F RELDF PRELDF RELDX STPPAR D*STEP NPRELDF
## 0 1 -4.227e+03
## 1 8 -4.227e+03 9.23e-06 1.81e-05 4.1e-05 4.5e+09 4.1e-06 4.05e+04
## 2 16 -4.231e+03 7.81e-04 1.35e-03 3.5e-01 2.0e+00 5.4e-02 1.33e+00
## 3 18 -4.233e+03 5.24e-04 6.88e-04 6.7e-01 1.9e+00 2.2e-01 2.15e-02
## 4 20 -4.233e+03 1.30e-04 1.75e-04 1.6e-01 8.9e-01 1.0e-01 2.79e-04
## 5 21 -4.234e+03 5.15e-05 1.14e-04 1.2e-01 9.7e-01 1.0e-01 2.11e-04
## 6 22 -4.234e+03 4.41e-05 2.78e-05 2.7e-02 0.0e+00 2.8e-02 2.78e-05
## 7 24 -4.234e+03 7.27e-05 7.72e-05 1.2e-01 2.3e-01 1.4e-01 7.97e-05
## 8 25 -4.234e+03 2.66e-05 1.88e-05 8.4e-02 0.0e+00 1.2e-01 1.88e-05
## 9 27 -4.234e+03 2.19e-05 1.84e-05 1.5e-02 2.0e+00 2.4e-02 6.19e-03
## 10 29 -4.235e+03 1.05e-04 6.54e-05 2.9e-02 2.0e+00 4.7e-02 1.75e+00
## 11 31 -4.237e+03 4.18e-04 3.42e-04 5.4e-02 2.0e+00 9.4e-02 4.39e+02
## 12 38 -4.237e+03 4.82e-06 9.44e-06 1.3e-07 5.1e+01 2.4e-07 5.37e-04
## 13 39 -4.237e+03 8.25e-09 2.17e-08 1.3e-07 2.0e+00 2.4e-07 1.16e-03
## 14 48 -4.238e+03 3.41e-04 6.41e-04 8.4e-03 8.4e-01 1.6e-02 1.16e-03
## 15 50 -4.239e+03 1.49e-04 1.51e-04 7.8e-03 1.5e+00 1.6e-02 2.72e-03
## 16 51 -4.240e+03 2.44e-04 3.34e-04 1.5e-02 1.8e+00 3.1e-02 3.90e-03
## 17 53 -4.240e+03 1.36e-05 1.54e-04 1.5e-03 2.0e+00 3.1e-03 2.50e-03
## 18 54 -4.240e+03 8.33e-05 9.68e-05 7.6e-04 2.0e+00 1.6e-03 4.03e-04
## 19 55 -4.240e+03 8.56e-06 1.01e-05 7.5e-04 2.0e+00 1.6e-03 6.28e-04
## 20 56 -4.240e+03 3.02e-06 5.67e-06 7.5e-04 1.9e+00 1.6e-03 1.42e-04
## 21 65 -4.240e+03 4.96e-07 9.56e-07 1.4e-08 5.4e+00 2.7e-08 1.10e-05
## 22 81 -4.240e+03 -8.79e-15 4.76e-14 1.0e-14 8.6e+04 2.0e-14 3.83e-06
##
## ***** FALSE CONVERGENCE *****
##
## FUNCTION -4.240102e+03 RELDX 1.036e-14
## FUNC. EVALS 81 GRAD. EVALS 22
## PRELDF 4.762e-14 NPRELDF 3.829e-06
##
## I FINAL X(I) D(I) G(I)
##
## 1 4.369506e-06 1.000e+00 1.003e+04
## 2 1.681183e-02 1.000e+00 -4.154e+01
## 3 9.723224e-01 1.000e+00 -2.010e+01fit$coe[3]## b1
## 0.9723224#so take the value of lambda as 0.97
#S&P500 index price
lp.GSPC2 = log(GSPC2[,5])
r11 = diff(lp.GSPC2)
r11 = r11[2:length(r11)]
fit = garch(r11)##
## ***** ESTIMATION WITH ANALYTICAL GRADIENT *****
##
##
## I INITIAL X(I) D(I)
##
## 1 5.092667e-05 1.000e+00
## 2 5.000000e-02 1.000e+00
## 3 5.000000e-02 1.000e+00
##
## IT NF F RELDF PRELDF RELDX STPPAR D*STEP NPRELDF
## 0 1 -5.452e+03
## 1 8 -5.454e+03 2.90e-04 4.93e-04 3.4e-05 2.4e+11 3.4e-06 5.84e+07
## 2 9 -5.454e+03 4.82e-06 5.72e-06 3.3e-05 2.0e+00 3.4e-06 1.60e+01
## 3 18 -5.477e+03 4.27e-03 6.93e-03 4.8e-01 2.0e+00 9.1e-02 1.58e+01
## 4 21 -5.512e+03 6.29e-03 5.75e-03 7.5e-01 1.9e+00 3.5e-01 7.23e-01
## 5 23 -5.519e+03 1.27e-03 2.37e-03 7.9e-02 2.0e+00 6.9e-02 2.15e+01
## 6 25 -5.536e+03 3.00e-03 4.87e-03 1.8e-01 2.0e+00 2.1e-01 2.38e+01
## 7 36 -5.541e+03 1.03e-03 2.84e-03 1.5e-06 4.1e+00 2.0e-06 4.56e-01
## 8 37 -5.542e+03 8.71e-05 6.60e-05 1.4e-06 2.0e+00 2.0e-06 2.94e-02
## 9 38 -5.542e+03 8.38e-06 9.86e-06 1.5e-06 2.0e+00 2.0e-06 1.81e-03
## 10 39 -5.542e+03 1.60e-07 1.14e-07 1.5e-06 2.0e+00 2.0e-06 5.15e-04
## 11 49 -5.542e+03 -1.17e-14 1.93e-15 3.4e-15 9.0e-01 4.6e-15 -2.62e-04
##
## ***** FALSE CONVERGENCE *****
##
## FUNCTION -5.541883e+03 RELDX 3.376e-15
## FUNC. EVALS 49 GRAD. EVALS 11
## PRELDF 1.931e-15 NPRELDF -2.617e-04
##
## I FINAL X(I) D(I) G(I)
##
## 1 6.141646e-06 1.000e+00 2.325e+03
## 2 2.200237e-01 1.000e+00 3.158e+01
## 3 6.817253e-01 1.000e+00 -5.507e+00fit$coe[3]## b1
## 0.6817253#so take the value of lambda as 0.68Define function of ordinary least squares based algorithm
- Define functions: Back to the main body
#model for exponential smoothing
vol.exp.sm <- function(x, lambda) {
sigma2 <- x^2
n <- length(x)
for (i in 2:n)
sigma2[i] <- sigma2[i-1] * lambda + x[i-1]^2 * (1-lambda)
sigma <- sqrt(sigma2)
resid <- x/sigma
resid[is.na(resid)] <- 0
sq.resid <- resid^2
list(sigma2=sigma2, sigma=sigma, resid = resid, sq.resid = sq.resid)
}
#model that computes the volatility for each covariate and then computes the acfs of the squared residuals after removing the volatility, and adds up
first.acf.squares.train <- function(x, lambda) {
d <- dim(x$x.train)
ss <- 0
x.train.dev <- x$x.train
y.train.dev <- x$y.train
x.valid.dev <- x$x.valid
y.valid.dev <- x$y.valid
x.test.dev <- x$x.test
y.test.dev <- x$y.test
for (i in 1:(d[2])) {
v <- vol.exp.sm(x$x.train[,i], lambda[i])
ss <- ss + abs(acf(v$sq.resid, plot=FALSE)$acf[2])
x.train.dev[,i] <- v$resid
v <- vol.exp.sm(x$x.valid[,i], lambda[i])
x.valid.dev[,i] <- v$resid
v <- vol.exp.sm(x$x.test[,i], lambda[i])
x.test.dev[,i] <- v$resid
}
v <- vol.exp.sm(x$y.train, lambda[i])
ss <- ss + abs(acf(v$sq.resid, plot=FALSE)$acf[2])
y.train.dev <- v$resid
v <- vol.exp.sm(x$y.valid, lambda[i])
y.valid.dev <- v$resid
v <- vol.exp.sm(x$y.test, lambda[i])
y.test.dev <- v$resid
list(ss=ss, y.train.dev=y.train.dev, x.train.dev=x.train.dev, y.valid.dev=y.valid.dev, x.valid.dev=x.valid.dev, y.test.dev=y.test.dev, x.test.dev=x.test.dev)
}
#estimation of beta in y = a + x beta + epsilon
thresh.reg <- function(x, y, th, x.pred = NULL) {
d <- dim(x)
ind <- (abs(cor(x, y)) > th)
n <- sum(ind)
new.x <- matrix(c(rep(1, d[1]), x[,ind]), d[1], n+1)
gram = t(new.x) %*% new.x
beta <- solve(gram) %*% t(new.x) %*% matrix(y, d[1], 1)
ind.ex <- c(1, as.numeric(ind))
ind.ex[ind.ex == 1] <- beta
condnum = max(svd(gram)$d)/min(svd(gram)$d)
pr <- 0
if (!is.null(x.pred)) pr <- sum(ind.ex * c(1, x.pred))
list(beta = ind.ex, pr=pr, condnum = condnum)
}
sharpe.curves <- function(x, lambda, th, warmup, reg.function = thresh.reg, win = seq(from = 10, to = warmup, by = 10)) {
w <- length(win)
train.curve <- valid.curve <- test.curve <- rep(0, w)
n <- length(x$y.train)
i=1
rreg <- rolling.thresh.reg(x, lambda, th, win[i], warmup, reg.function)
rreg.valid <- rolling.thresh.reg.valid(x, lambda, th, win[i], warmup, reg.function)
rreg.test <- rolling.thresh.reg.test(x, lambda, th, win[i], warmup, reg.function)
condnum = matrix(0,w,length(rreg$condnum))
condnum.valid = matrix(0,w,length(rreg.valid$condnum))
condnum.test = matrix(0,w,length(rreg.test$condnum))
train.curve[i] <- rreg$err
valid.curve[i] <- rreg.valid$err
test.curve[i] <- rreg.test$err
condnum[i,] <- rreg$condnum
condnum.valid[i,] <- rreg.valid$condnum
condnum.test[i,] <- rreg.test$condnum
for (i in 2:w) {
rreg <- rolling.thresh.reg(x, lambda, th, win[i], warmup, reg.function)
rreg.valid <- rolling.thresh.reg.valid(x, lambda, th, win[i], warmup, reg.function)
rreg.test <- rolling.thresh.reg.test(x, lambda, th, win[i], warmup, reg.function)
train.curve[i] <- rreg$err
valid.curve[i] <- rreg.valid$err
test.curve[i] <- rreg.test$err
condnum[i,] <- rreg$condnum
condnum.valid[i,] <- rreg.valid$condnum
condnum.test[i,] <- rreg.test$condnum
}
list(train.curve = train.curve, valid.curve = valid.curve, test.curve = test.curve, condnum=condnum, condnum.valid = condnum.valid, condnum.test = condnum.test)
}
rolling.thresh.reg <- function(x, lambda, th, win, warmup, reg.function = thresh.reg) {
xx <- first.acf.squares.train(x, lambda)
n <- length(xx$y.train.dev)
err <- 0
condnum <- predi <- truth <- rep(0, n-warmup+1)
for (i in warmup:n) {
y <- xx$y.train.dev[(i-win):(i-1)]
xxx <- xx$x.train.dev[(i-win):(i-1),]
zz <- reg.function(xxx, y, th, xx$x.train.dev[i,])
predi[i-warmup+1] <- zz$pr
condnum[i-warmup+1] <- zz$condnum
truth[i-warmup+1] <- xx$y.train.dev[i]
}
predi = sign(predi) *1* (abs(predi) > 0)
ret <- predi * truth
err <- sqrt(250) * mean(ret) / sqrt(var(ret))
list(err=err, predi=predi, truth=truth, condnum=condnum)
}
rolling.thresh.reg.valid <- function(x, lambda, th, win, warmup, reg.function = thresh.reg) {
xx <- first.acf.squares.train(x, lambda)
n <- length(xx$y.valid.dev)
err <- 0
condnum <- predi <- truth <- rep(0, n-warmup+1)
for (i in warmup:n) {
y <- xx$y.valid.dev[(i-win):(i-1)]
xxx <- xx$x.valid.dev[(i-win):(i-1),]
zz <- reg.function(xxx, y, th, xx$x.valid.dev[i,])
predi[i-warmup+1] <- zz$pr
condnum[i-warmup+1] <- zz$condnum
truth[i-warmup+1] <- xx$y.valid.dev[i]
}
predi = sign(predi) *1* (abs(predi) > 0)
ret <- predi * truth
err <- sqrt(250) * mean(ret) / sqrt(var(ret))
list(err=err, predi=predi, truth=truth, condnum=condnum)
}
rolling.thresh.reg.test <- function(x, lambda, th, win, warmup, reg.function = thresh.reg) {
xx <- first.acf.squares.train(x, lambda)
n <- length(xx$y.test.dev)
err <- 0
condnum <- predi <- truth <- rep(0, n-warmup+1)
for (i in warmup:n) {
y <- xx$y.test.dev[(i-win):(i-1)]
xxx <- xx$x.test.dev[(i-win):(i-1),]
zz <- reg.function(xxx, y, th, xx$x.test.dev[i,])
predi[i-warmup+1] <- zz$pr
condnum[i-warmup+1] <- zz$condnum
truth[i-warmup+1] <- xx$y.test.dev[i]
}
predi = sign(predi) *1* (abs(predi) > 0)
ret <- predi * truth
err <- sqrt(250) * mean(ret) / sqrt(var(ret))
list(err=err, predi=predi, truth=truth, condnum=condnum)
}Detailed explanation of concpet of ordinary least squares based algorithm
- Concept for the algorithm: Back to the main body
First, define “local” data using a tuning parameter D, representing the number of past data points using for estimating the model at current time t.
To simplify further the notations, the local model at time t: \[ Y_{s+1} = \alpha_0+\sum_{k\ge0}^{q+10}\alpha_kZ^k_s+e_{s+1}, \ s = t-D, ..., t-1 \] It uses a warmup time t0, a window length D, and the 10+q return series as input. The model above could be rewritten as Y = Zα +e, where \[ \mathbf{Y} = (Y_{t-D+1}, ...,Y_t)^T,\\ \mathbf{Z} = (\mathbf{1}_D, \mathbf{Z}^1, ..., \mathbf{Z}^{q+10})\\ \text{with} ~ \mathbf{Z}^k = (Z^{k}_{t-D}, ..., Z^{k}_{t-1})^T, ~\alpha = (\alpha_0, \alpha_1, ...,\alpha_{q+10})^T ~\text{and}~ \mathbf{e} = (e_{t-D+1},..., e_t)^T \] The predicted return at time t+1 is then \[ \hat Y_{t+1} = \hat\alpha^T\mathbf{z}^t, ~ \text{where}~ \mathbf{z}_t = (1, Z^1_t,...,Z^{q+10}_t)^T \] The model also considers the next day’s S&P500 return. If the return is positive, 1 unit of money should be invested into S&P500, and it would invest -1 unit of money if the prediction is negative.
It should be mentioned that the model is using the short_selling strategy, which is that one unit of money is invested into the market at the start of the day and only read the predicted return before closing. If the next day’s return is predicted to be negative, the money is no longer kept in the market but back to hand. If the next day’s return is predicted to be positive, the money is saved in the trading account for the next day.
Hence I made additional notes in three rolling.thresh.reg functions.
I consider all negative values as -1 and all positive values as 1. So
the code predi = sign(predi) *1* (abs(predi) > 0) is
added before defining the actual return and the calculation of the
Sharpe ratio.
Take the predicted return as the new position (ignore all transaction costs to achieve this). I simplify the procedure by not multiplying back the forecasted volatility effect at time t+1. With this, the actual return at time t+1 is: (with the adjusted predicted return)
\[ R_{t+1} = \hat{Y}_{t+1}Y_{t+1} \] And then the whole process is started again but with t replaced by t +1. Stop when t reaches the end of the training set and repeat it for a grid of values of D.
This tuning parameter D is to be determined using the validation set, and I use Sharpe ratio to see which D is doing the best job under the training set.
Sharpe ratio is defined by: \[ \text{Sharpe ratio} = \sqrt{250}\times\frac{\text{Average daily return}}{\text{SD of daily returns}} \] As in former process, the tuning parameter D is directly related to the fact that the model is time dependent. The same model trained with a particular D in the training set cannot be used in the validation set, but has to be estimated from the validation data again.
So the validation procedure is simply to repeat the steps for the validation set, producing different Sharpe ratios for different values of D, then compare them to find which D is the best from the training and validation sets.
Define function of Principal component regression based algorithm
- Set GARCH(1,1) model for each variables: Back to the main body
#new function
rolling.thresh.reg <- function(x, lambda, th, win, warmup, reg.function = thresh.reg, PCR_analysis = FALSE) {
xx <- first.acf.squares.train(x, lambda)
if(PCR_analysis){
train.pcr = prcomp(x$x.train)
x$x.train = train.pcr$x[,1:2]
}
n <- length(xx$y.train.dev)
err <- 0
condnum <- predi <- truth <- rep(0, n-warmup+1)
for (i in warmup:n) {
y <- xx$y.train.dev[(i-win):(i-1)]
xxx <- xx$x.train.dev[(i-win):(i-1),]
zz <- reg.function(xxx, y, th, xx$x.train.dev[i,])
predi[i-warmup+1] <- zz$pr
condnum[i-warmup+1] <- zz$condnum
truth[i-warmup+1] <- xx$y.train.dev[i]
}
predi = sign(predi) *1* (abs(predi) > 0)
ret <- predi * truth
err <- sqrt(250) * mean(ret) / sqrt(var(ret))
list(err=err, predi=predi, truth=truth, condnum=condnum)
}
rolling.thresh.reg.valid <- function(x, lambda, th, win, warmup, reg.function = thresh.reg, PCR_analysis = FALSE) {
xx <- first.acf.squares.train(x, lambda)
if(PCR_analysis){
valid.pcr = prcomp(x$x.valid)
x$x.valid = valid.pcr$x[,1:2]
}
n <- length(xx$y.valid.dev)
err <- 0
condnum <- predi <- truth <- rep(0, n-warmup+1)
for (i in warmup:n) {
y <- xx$y.valid.dev[(i-win):(i-1)]
xxx <- xx$x.valid.dev[(i-win):(i-1),]
zz <- reg.function(xxx, y, th, xx$x.valid.dev[i,])
predi[i-warmup+1] <- zz$pr
condnum[i-warmup+1] <- zz$condnum
truth[i-warmup+1] <- xx$y.valid.dev[i]
}
predi = sign(predi) *1* (abs(predi) > 0)
ret <- predi * truth
err <- sqrt(250) * mean(ret) / sqrt(var(ret))
list(err=err, predi=predi, truth=truth, condnum=condnum)
}
rolling.thresh.reg.test <- function(x, lambda, th, win, warmup, reg.function = thresh.reg, PCR_analysis = FALSE) {
xx <- first.acf.squares.train(x, lambda)
if(PCR_analysis){
test.pcr = prcomp(x$x.test)
x$x.test = test.pcr$x[,1:2]
}
n <- length(xx$y.test.dev)
err <- 0
condnum <- predi <- truth <- rep(0, n-warmup+1)
for (i in warmup:n) {
y <- xx$y.test.dev[(i-win):(i-1)]
xxx <- xx$x.test.dev[(i-win):(i-1),]
zz <- reg.function(xxx, y, th, xx$x.test.dev[i,])
predi[i-warmup+1] <- zz$pr
condnum[i-warmup+1] <- zz$condnum
truth[i-warmup+1] <- xx$y.test.dev[i]
}
predi = sign(predi) *1* (abs(predi) > 0)
ret <- predi * truth
err <- sqrt(250) * mean(ret) / sqrt(var(ret))
list(err=err, predi=predi, truth=truth, condnum=condnum)
}Importing xts objects
- Importing data - 1: Back to the main body
getSymbols(c("AAPL","MSFT","AMZN","GOOG","AVGO","TSLA","NVDA","PEP", "META","COST","^GSPC"), src = "yahoo", from = as.Date("2013-01-01"), to = as.Date("2017-12-01"))## [1] "AAPL" "MSFT" "AMZN" "GOOG" "AVGO" "TSLA" "NVDA" "PEP" "META"
## [10] "COST" "^GSPC"AAPL2 <- AAPL
MSFT2 <- MSFT
AMZN2 <- AMZN
GOOG2 <- GOOG
AVGO2 <- AVGO
TSLA2 <- TSLA
NVDA2 <- NVDA
PEP2 <- PEP
META2 <- META
COST2 <- COST
GSPC2 <- GSPC
getSymbols(c("AAPL","MSFT","AMZN","GOOG","AVGO","TSLA","NVDA","PEP", "META","COST","^GSPC"), src = "yahoo", from = as.Date("2018-01-01"), to = as.Date("2022-12-01"))## [1] "AAPL" "MSFT" "AMZN" "GOOG" "AVGO" "TSLA" "NVDA" "PEP" "META"
## [10] "COST" "^GSPC"str(AAPL)## An 'xts' object on 2018-01-02/2022-11-30 containing:
## Data: num [1:1238, 1:6] 42.5 43.1 43.1 43.4 43.6 ...
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : chr [1:6] "AAPL.Open" "AAPL.High" "AAPL.Low" "AAPL.Close" ...
## Indexed by objects of class: [Date] TZ: UTC
## xts Attributes:
## List of 2
## $ src : chr "yahoo"
## $ updated: POSIXct[1:1], format: "2023-03-06 00:38:31"str(AAPL2)## An 'xts' object on 2013-01-02/2017-11-30 containing:
## Data: num [1:1239, 1:6] 19.8 19.6 19.2 18.6 18.9 ...
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : chr [1:6] "AAPL.Open" "AAPL.High" "AAPL.Low" "AAPL.Close" ...
## Indexed by objects of class: [Date] TZ: UTC
## xts Attributes:
## List of 2
## $ src : chr "yahoo"
## $ updated: POSIXct[1:1], format: "2023-03-06 00:38:22"Detailed explanation of concept for Value-at-Risk - 1
- Using ARIMA model to estimate Value-at-Risk: Back to the main body
Non-stationary processes have means, variances and covariances that change over time. Using non-stationary time series data leads to unreliable forecasting. A stationary process is mean-reverting, i.e, it fluctuates around a constant mean with constant variance.
In this case, stationarity refers to weak stationarity where the stationary time series satisfies three conditions \[ E(y_t) = \mu\\ E(y^2_t) = \sigma^2 < \infty\\ Cov(y_1, y_{1+s}) = Cov(y_2, y_{2+s}) = \cdots = Cov(y_T, y_{T+s}) = \gamma_s \] In order to resolve this issue, the 1-order differencing can be described as \[ \Delta y_t = y_t - y_{t-1} \] For the stationarity transformation, we prefer to calculate the simple daily returns, expressed as follows \[ r_t = \frac{price_t - price_{t-1}}{price_{t-1}} \] The Box-Jenkins approach is commonly used for time series analysis. It involves the application of ARIMA models to identify the optimal time series model that represents the underlying stochastic process.
Before implementing this methodology, it is necessary to ensure that the time series is stationary. In this section, I utilized stock returns that had previously been evaluated for stationarity. The identification of the ARIMA model order can be determined using the Autocorrelation Function (ACF), Partial Autocorrelation Function (PACF), or the Akaike Information Criterion (AICc). The AICc estimates the quality of each model relative to one another.
\[
AIC = \ln\left(\frac{\sum_{i=1}^{T} u_i^2}{T}\right) + \frac{2k}{T}
\] where \sum_{i=1}^{T} u_i^2 is sum of Squared
Residuals; T is number of observations, and k
is number of model parameters (p+q+1).
The addition of extra lag parameters to a model can result in a decrease in the Sum of Squared Residuals. However, this approach may also lead to overfitting. The Akaike Information Criterion (AIC) addresses the risk of both overfitting and underfitting, so the model with the lowest AIC value is considered the best fit.
Detailed explanation of concept for Value-at-Risk - 2
- Using GARCH(1,1) to estimate Value-at-Risk: Back to the main body
Value at risk, or VaR is a statistic that quantifies the extent of possible financial losses within a firm, portfolio, or position over a specific time frame. It is an estimate of the maximum amount to lose with a certain probability during a particular holding period.
VaR is mathematically defined as \[ \text{VaR}(X) = \mu + \sigma{N}^{-1}(X) \] where is the mean of daily returns, is the (usually sample) standard deviation representing the unconditional volatility and {N}^{-1} is the inverse of the cumulative normal distribution at a given confidence level X.
Since the mean is approximately zero, the formula below is used to calculate the volatility \[ \sigma^2_{n} = \frac{1}{n}\sum_{i=1}^n R_{i}^2 \] where n is the number of observations. This formula ignores “volatility clustering”. By doing so, risk is underestimated/overestimated during periods with high/low volatility periods.
In order to incorporate the conditional volatility into the VaR formula, I use a GARCH(1,1), which is defined as \[ \sigma^2_t = \omega + \alpha R_{t-1}^2 + \beta \sigma^2_{t-1} \]