#Run installed library library(TSA) library(tseries) library(rugarch) library(quantmod) library(zoo) library(xts) library(MASS) library(forecast) library(dplyr)

#Get data from Yahoo #Fix the date period start <- as.Date(“2017-10-01”) end <- as.Date(“2020-10-31”) RDS <- getSymbols(“RDS-A”, src = “yahoo”, from = start, to = end, warnings = FALSE, auto.assign = TRUE)

#Data Summary head(RDS-A) summary(RDS-A)

#Assign names Price <- quantmod::Ad(RDS-A) head(Price) summary(Price) Date <- index(RDS-A) head(Date)

#Plot for the original dataset plot(Price, ylim = c(10,80), main = “Royal Dutch Shell Stock Price”)

#Plot for log price logprice <- log(Price) plot(logprice, ylim = c(2.5, 4.5), main = “Royal Dutch Shell Stock Log Price”) acf(logprice,50) pacf(logprice,50) adf.test(logprice, alternative = “stationary”) summary(Price)

#Apply first differencing return <- diff(logprice) plot(return, main = “RDS Daily Return”) summary(return) acf(return, lag.max = 50, na.action = na.pass, main = “ACF of Daily Return”) pacf(return, lag.max = 50, na.action = na.pass, main = “PACF of Daily Return”) kpss.test(return)

#Linear Models #Parameter Estimation + Diagnostic Checking Model1 <- Arima(return, order = c(2,0,2), method = “ML”) summary(Model1)

Model2 <- Arima(logprice, order = c(2,1,2), method = “ML”) summary(Model2)

Model3 <- Arima(logprice, order = c(6,1,6), method = “ML”) summary(Model3)

Model4 <- Arima(logprice, order = c(7,1,7), method = “ML”) summary(Model4)

Model5 <- Arima(logprice, order = c(8,1,8), method = “ML”) summary(Model5)

#Non-linear Model #Fit the model with white noise and find residual white.noise <- Arima(return, order= c(0,0,0)) summary(white.noise) checkresiduals(white.noise)

residual <- white.noise$residuals s.residual <- residual^2

par(mfcol = c(3,1)) acf(residual, max.lag = 30, na.action = na.pass, main = “ACF of Residuals”) pacf(residual, max.lag = 30, na.action = na.pass, main = “PACF of Residuals”) plot(residual, main = “Residuals”)

par(mfcol = c(3,1)) acf(s.residual, max.lag = 30, na.action = na.pass, main = “ACF of Squared Residuals”) pacf(s.residual, max.lag = 30, na.action = na.pass, main = “PACF of Squared Residuals”) plot(s.residual, main = “Squared Residual”)

#ARCH Modelling #Model A : GARCH(1,0) ARMA(0,0) garch.1000 <- ugarchspec(variance.model = list(model = “sGARCH”, garchOrder = c(1,0)), mean.model = list(armaOrder = c(0,0)), distribution.model = “std”) ModelA <- ugarchfit(spec = garch.1000, data = logprice) ModelA

#Model B : GARCH(1,0) ARMA(1,0) garch.1010 <- ugarchspec(variance.model = list(model = “sGARCH”, garchOrder = c(1,0)), mean.model = list(armaOrder = c(1,0)), distribution.model = “std”) ModelB <- ugarchfit(spec = garch.1010, data = logprice) ModelB

#Model C : GARCH(1,0) ARMA(2,2) garch.1022 <- ugarchspec(variance.model = list(model = “sGARCH”, garchOrder = c(1,0)), mean.model = list(armaOrder = c(2,2)), distribution.model = “std”) ModelC <- ugarchfit(spec = garch.1022, data = logprice) ModelC

#GARCH Modelling #Model D : GARCH(1,1) ARMA(0,0) garch.1100 <- ugarchspec(variance.model = list(model = “sGARCH”, garchOrder = c(1,1)), mean.model = list(armaOrder = c(0,0)), distribution.model = “std”) ModelD <- ugarchfit(spec = garch.1100, data = logprice) ModelD

#Model E : GARCH(1,1) ARMA(1,0) garch.1110 <- ugarchspec(variance.model = list(model = “sGARCH”, garchOrder = c(1,1)), mean.model = list(armaOrder = c(1,0)), distribution.model = “std”) ModelE <- ugarchfit(spec = garch.1110, data = logprice) ModelE

#Model F : GARCH(1,1) ARMA(2,2) garch.1122 <- ugarchspec(variance.model = list(model = “sGARCH”, garchOrder = c(1,1)), mean.model = list(armaOrder = c(2,2)), distribution.model = “std”) ModelF <- ugarchfit(spec = garch.1122, data = logprice) ModelF

#Forecasting on Linear Model - Model 2 forecast10.model2 <- forecast(Model2,10, level = 95) #10-step forecast autoplot(forecast10.model2, series = “Prediction”, xlim = c(750,800)) forecast10.model2

#Forecasting on non-Linear Model forecast.modelE <- ugarchboot(ModelE, n.ahead = 10, method = c(“Partial”, “Full”)[1]) plot(forecast.modelE, which = 2) forecast.modelE

forecast.modelF <- ugarchboot(ModelF, n.ahead = 10, method = c(“Partial”, “Full”)[1]) plot(forecast.modelF, which = 2) forecast.modelF