Forecasting Assignment 1
Student Details
- Shipren Jayadev (s3744421)
Data
The data is obtained from the ASX historical data which covers monthly averages of ASX All Ordinaries (Ords) Price Index, Gold price (AUD), Crude Oil (Brent, USD/bbl) and Copper (USD/tonne) starting from January 2004.
Introduction
The purpose of this report is to analyse the dataset and be able to determine the following:
- The existence of nonstationarity in the dataset.
- The impact of the components of a time series data on the given dataset.
- The most accurate and suitable distributed lag model for the ASX price index.
Load packages
library(readr)
library(tseries)
library(forecast)
library(TSA)
library(x12)
library(dLagM)
library(car)Load data
asx <- read_csv("~/RMIT Classes/Semester 3/Forecasting/Assignment 1/ASX_data.csv")Processing and visual description of data
#Check the classes of the data
head(asx)class(asx)## [1] "spec_tbl_df" "tbl_df" "tbl" "data.frame"#Convert data to time series
price.ts <- ts(asx$`ASX price`,start= c(2004, 1), frequency = 12)
gold.ts <- ts(asx$`Gold price`, start= c(2004, 1), frequency = 12)
oil.ts <- ts(asx$`Crude Oil (Brent)_USD/bbl`, start= c(2004, 1), frequency = 12)
copper.ts <- ts(asx$`Copper_USD/tonne`, start= c(2004, 1), frequency = 12)
#Plot the new time series variables
plot(price.ts)
plot(gold.ts)
plot(oil.ts)
plot(copper.ts)
Checking for Non-stationarity
#Plot the ACF and PACF for each variable
par(mfrow=c(1,2))
acf(price.ts)
pacf(price.ts)
acf(gold.ts)
pacf(gold.ts)
acf(oil.ts)
pacf(oil.ts)
acf(copper.ts)
pacf(copper.ts)
#Estimating lambda for each variable
lambda.price <- BoxCox.lambda(price.ts) %>% print()## [1] 1.999924lambda.gold <- BoxCox.lambda(gold.ts) %>% print()## [1] 0.976695lambda.oil <- BoxCox.lambda(oil.ts) %>% print()## [1] -0.8304136lambda.copper <- BoxCox.lambda(copper.ts) %>% print()## [1] 0.9336783#No transformation needed for copper and gold variables since lambda is close to 1
#Applying box-cox transformation to oil and price
bc.price <- BoxCox(price.ts, lambda = lambda.price)
bc.oil <- BoxCox(oil.ts, lambda = lambda.oil)
#Conducting ADF test for transformed price and oil variables
k.price <- ar(bc.price)$order
adf.test(bc.price, k = k.price)##
## Augmented Dickey-Fuller Test
##
## data: bc.price
## Dickey-Fuller = -2.1072, Lag order = 1, p-value = 0.5316
## alternative hypothesis: stationaryk.oil <- ar(bc.oil)$order
adf.test(bc.oil, k = k.oil)##
## Augmented Dickey-Fuller Test
##
## data: bc.oil
## Dickey-Fuller = -1.7787, Lag order = 1, p-value = 0.6686
## alternative hypothesis: stationary#Results still indicate the variables are nonstationary
#Applying differencing transformation
par(mfrow = c(1,1))
price.diff <- diff(bc.price)
plot(price.diff, main = 'First Difference Price Time Series Plot')
gold.diff <- diff(gold.ts)
plot(gold.diff, main = 'First Difference Gold Time Series Plot')
oil.diff <- diff(bc.oil)
plot(price.diff, main = 'First Difference Oil Time Series Plot')
copper.diff <- diff(copper.ts)
plot(oil.diff, main = 'First Difference Copper Time Series Plot')
#All variables appear to be stationary
#Apply ADF test to ensure time series is stationary
k2.price <- ar(price.diff)$order
adf.test(price.diff, k = k2.price)## Warning in adf.test(price.diff, k = k2.price): p-value smaller than printed p-
## value##
## Augmented Dickey-Fuller Test
##
## data: price.diff
## Dickey-Fuller = -11.692, Lag order = 0, p-value = 0.01
## alternative hypothesis: stationaryk2.gold <- ar(gold.diff)$order
adf.test(gold.diff, k = k2.gold)## Warning in adf.test(gold.diff, k = k2.gold): p-value smaller than printed p-
## value##
## Augmented Dickey-Fuller Test
##
## data: gold.diff
## Dickey-Fuller = -7.3956, Lag order = 2, p-value = 0.01
## alternative hypothesis: stationaryk2.oil <- ar(oil.diff)$order
adf.test(oil.diff, k = k2.oil)## Warning in adf.test(oil.diff, k = k2.oil): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: oil.diff
## Dickey-Fuller = -7.9276, Lag order = 1, p-value = 0.01
## alternative hypothesis: stationaryk2.copper <- ar(copper.diff)$order
adf.test(copper.diff, k = k2.copper)## Warning in adf.test(copper.diff, k = k2.copper): p-value smaller than printed p-
## value##
## Augmented Dickey-Fuller Test
##
## data: copper.diff
## Dickey-Fuller = -7.3527, Lag order = 1, p-value = 0.01
## alternative hypothesis: stationary#ADF test now indicates variables are stationary##Existence of Nonstationarity
When first conducting the analysis, the ACF and PACF graphs indicate the time series are nonstationary and there is no indication of seasonality but there is an observable downward trend.
The lambda for each time series is obtained. Time series with lambda not close to 1 undergo a boxcox transformation and then an ADF test to determine if it has become stationary.
Price and oil had lambda values not close to 1 so underwent a boxcox transformation and the ADF test determined the variables were still nonstationary.
All 4 variables underwent differencing and an ADF test was run on all 4 variables which determined they were now all stationary.
Decomposition of Time Series Data
#Applying x12 decomposition methods to observe the effects of trend and seasonality
price.window <- window(bc.price, start = c(2004,1))
price.decom <- x12(price.window)
plot(price.decom, sa = TRUE, trend = TRUE, forecast = TRUE)
plotSeasFac(price.decom, main = 'Price Time Series Seasonal Factors')
gold.window <- window(gold.ts, start = c(2004,1))
gold.decom <- x12(gold.window)
plot(gold.decom, sa = TRUE, trend = TRUE, forecast = TRUE)
plotSeasFac(gold.decom, main = 'Gold Time Series Seasonal Factors')
oil.window <- window(bc.oil, start = c(2004,1))
oil.decom <- x12(oil.window)
plot(oil.decom, sa = TRUE, trend = TRUE, forecast = TRUE)
plotSeasFac(oil.decom, main = 'Oil Time Series Seasonal Factors')
copper.window <- window(copper.ts, start = c(2004,1))
copper.decom <- x12(copper.window)
plot(copper.decom, sa = TRUE, trend = TRUE, forecast = TRUE)
plotSeasFac(copper.decom, main = 'Copper Time Series Seasonal Factors')
Impact of Components on the Data
Price: There is very little deviation between the original and seasonally adjusted line, implying no seasonality as initially stated. There is an intervention point in the year 2008. There are not a lot of extreme values but there are large deviations from the mean in the months of February, July and September.
Gold: There is a sudden upward spike in the year 2006 indicating an intervention point. Seasonality has very little effect on the series as it does not deviate much from the original line. There are a large deviations from the mean in the months of March to August, November and December.
Oil: There is a steep decline in 2004 followed by a sharp increase indicating an intervention point. The same can be said for the year 2008 and 2014 with large decreases. There are large deviations from the mean particularly in the months of June and October.
Copper: There is a sharp increase starting in 2006 implying intervention as well as a sharp decrease in 2008 and large increase startgin in 2009, also indicating intervention points.There are large fluctuations around the mean in all months except for May and June where there is little fluctuation around the mean.
Distributed lag model for ASX price index
#Plotting asx time series
asx.ts <- ts(asx, start = c(2004,1), frequency = 12)
asx.scaled = scale(asx.ts)
asxdata <- as.data.frame(asx)
plot(asx.scaled, plot.type="s",col = c("blue", "red", "green", "orange"), main = "Scaled ASX Variables")
legend("topleft",lty=1, text.width = 1.5, col=c("blue","red", "green", "orange"), c("Price", "Gold", "Oil", "Copper"))
cor(asx.ts)## ASX price Gold price Crude Oil (Brent)_USD/bbl
## ASX price 1.0000000 0.3431908 0.3290338
## Gold price 0.3431908 1.0000000 0.4366382
## Crude Oil (Brent)_USD/bbl 0.3290338 0.4366382 1.0000000
## Copper_USD/tonne 0.5617864 0.5364213 0.8664296
## Copper_USD/tonne
## ASX price 0.5617864
## Gold price 0.5364213
## Crude Oil (Brent)_USD/bbl 0.8664296
## Copper_USD/tonne 1.0000000DLM Fitting for Gold Price Variable
# Fitting a DLM model
for ( i in 1:6){
gold.model1.1 <- dlm(x = as.vector(asx$`Gold price`) , y = as.vector(asx$`ASX price`), q = i)
cat("q = ", i, "AIC = ", AIC(gold.model1.1$model), "BIC = ", BIC(gold.model1.1$model),"\n")
}## q = 1 AIC = 2613.609 BIC = 2625.91
## q = 2 AIC = 2596.292 BIC = 2611.637
## q = 3 AIC = 2579.215 BIC = 2597.59
## q = 4 AIC = 2562.296 BIC = 2583.69
## q = 5 AIC = 2544.887 BIC = 2569.286
## q = 6 AIC = 2527.575 BIC = 2554.966#Best DLM model with lag = 6
gold.model <- dlm( x = as.vector(asx$`Gold price`) , y = as.vector(asx$`ASX price`), q = 6)
summary(gold.model)##
## Call:
## lm(formula = model.formula, data = design)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1498.74 -617.99 -5.33 492.71 2084.31
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4250.19119 217.51364 19.540 <2e-16 ***
## x.t -0.19930 1.26861 -0.157 0.875
## x.1 -0.04654 1.90563 -0.024 0.981
## x.2 0.07110 1.94292 0.037 0.971
## x.3 -0.05378 1.95473 -0.028 0.978
## x.4 -0.11744 1.94739 -0.060 0.952
## x.5 0.27967 1.91454 0.146 0.884
## x.6 0.61039 1.26350 0.483 0.630
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 814.8 on 147 degrees of freedom
## Multiple R-squared: 0.07473, Adjusted R-squared: 0.03067
## F-statistic: 1.696 on 7 and 147 DF, p-value: 0.1142
##
## AIC and BIC values for the model:
## AIC BIC
## 1 2527.575 2554.966gold.VIF.model <- vif(gold.model$model)
gold.VIF.model## x.t x.1 x.2 x.3 x.4 x.5 x.6
## 56.34101 128.59349 134.96243 138.23712 138.55461 134.96320 59.24582#Fitting a ploynomial model
for (i in 1:6){
gold.model2.1 <- polyDlm(x = as.vector(asx$`Gold price`) , y = as.vector(asx$`ASX price`), q = i , k = i, show.beta = FALSE)
cat("q = ", i, "k = ", i, "AIC = ", AIC(gold.model2.1$model), "BIC = ", BIC(gold.model2.1$model),"\n")
}## q = 1 k = 1 AIC = 2613.609 BIC = 2625.91
## q = 2 k = 2 AIC = 2596.292 BIC = 2611.637
## q = 3 k = 3 AIC = 2579.215 BIC = 2597.59
## q = 4 k = 4 AIC = 2562.296 BIC = 2583.69
## q = 5 k = 5 AIC = 2544.887 BIC = 2569.286
## q = 6 k = 6 AIC = 2527.575 BIC = 2554.966#Best polyDLM model with q = 6 and i = 6
gold.model2 <- polyDlm(x = as.vector(asx$`Gold price`) , y = as.vector(asx$`ASX price`), q = 6, k = 6, show.beta = FALSE)
summary(gold.model2)##
## Call:
## "Y ~ (Intercept) + X.t"
##
## Residuals:
## Min 1Q Median 3Q Max
## -1498.74 -617.99 -5.33 492.71 2084.31
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.250e+03 2.175e+02 19.540 <2e-16 ***
## z.t0 -1.993e-01 1.269e+00 -0.157 0.875
## z.t1 -8.414e-03 3.749e+01 0.000 1.000
## z.t2 2.623e-01 7.118e+01 0.004 0.997
## z.t3 -9.368e-02 4.890e+01 -0.002 0.998
## z.t4 -1.579e-02 1.543e+01 -0.001 0.999
## z.t5 9.191e-03 2.269e+00 0.004 0.997
## z.t6 -8.435e-04 1.260e-01 -0.007 0.995
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 814.8 on 147 degrees of freedom
## Multiple R-squared: 0.07473, Adjusted R-squared: 0.03067
## F-statistic: 1.696 on 7 and 147 DF, p-value: 0.1142gold.VIF.model2 <- vif(gold.model2$model)
gold.VIF.model2## z.t0 z.t1 z.t2 z.t3 z.t4 z.t5
## 2.786620e+03 2.231986e+07 1.523293e+09 1.697529e+10 4.511506e+10 2.815051e+10
## z.t6
## 2.635256e+09#Applying koyck transformation
gold.model3 <- koyckDlm(x = as.vector(asx$`Gold price`) , y = as.vector(asx$`ASX price`))
summary(gold.model3)##
## Call:
## "Y ~ (Intercept) + Y.1 + X.t"
##
## Residuals:
## Min 1Q Median 3Q Max
## -682.19 -105.44 15.86 135.04 783.60
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.902e+02 8.958e+01 2.123 0.0353 *
## Y.1 9.635e-01 1.909e-02 50.469 <2e-16 ***
## X.t 2.595e-03 4.304e-02 0.060 0.9520
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 201.4 on 157 degrees of freedom
## Multiple R-Squared: 0.9488, Adjusted R-squared: 0.9481
## Wald test: 1454 on 2 and 157 DF, p-value: < 2.2e-16
##
## alpha beta phi
## Geometric coefficients: 5205.15 0.002595168 0.9634602gold.VIF.model3 <- vif(gold.model3$model)
gold.VIF.model3## Y.1 X.t
## 1.140648 1.140648#Fitting an ARDL model
for (i in 1:6){
for(j in 1:6){
gold.model4.1 = ardlDlm(x = as.vector(asx$`Gold price`) , y = as.vector(asx$`ASX price`), p = i , q = j )
cat("p = ", i, "q = ", j, "AIC = ", AIC(gold.model4.1$model), "BIC = ", BIC(gold.model4.1$model),"\n")
}
}## p = 1 q = 1 AIC = 2140.897 BIC = 2156.273
## p = 1 q = 2 AIC = 2128.524 BIC = 2146.938
## p = 1 q = 3 AIC = 2113.99 BIC = 2135.428
## p = 1 q = 4 AIC = 2102.754 BIC = 2127.204
## p = 1 q = 5 AIC = 2092.194 BIC = 2119.643
## p = 1 q = 6 AIC = 2081.15 BIC = 2111.584
## p = 2 q = 1 AIC = 2128.627 BIC = 2147.04
## p = 2 q = 2 AIC = 2130.523 BIC = 2152.005
## p = 2 q = 3 AIC = 2115.89 BIC = 2140.39
## p = 2 q = 4 AIC = 2104.694 BIC = 2132.2
## p = 2 q = 5 AIC = 2094.14 BIC = 2124.639
## p = 2 q = 6 AIC = 2083.058 BIC = 2116.535
## p = 3 q = 1 AIC = 2118.109 BIC = 2139.547
## p = 3 q = 2 AIC = 2120.027 BIC = 2144.528
## p = 3 q = 3 AIC = 2117.305 BIC = 2144.868
## p = 3 q = 4 AIC = 2105.731 BIC = 2136.293
## p = 3 q = 5 AIC = 2095.264 BIC = 2128.812
## p = 3 q = 6 AIC = 2084.237 BIC = 2120.758
## p = 4 q = 1 AIC = 2107.002 BIC = 2131.452
## p = 4 q = 2 AIC = 2108.914 BIC = 2136.42
## p = 4 q = 3 AIC = 2106.276 BIC = 2136.839
## p = 4 q = 4 AIC = 2107.456 BIC = 2141.074
## p = 4 q = 5 AIC = 2097.01 BIC = 2133.608
## p = 4 q = 6 AIC = 2085.895 BIC = 2125.46
## p = 5 q = 1 AIC = 2094.908 BIC = 2122.357
## p = 5 q = 2 AIC = 2096.86 BIC = 2127.359
## p = 5 q = 3 AIC = 2094.144 BIC = 2127.692
## p = 5 q = 4 AIC = 2095.425 BIC = 2132.023
## p = 5 q = 5 AIC = 2097.324 BIC = 2136.972
## p = 5 q = 6 AIC = 2086.36 BIC = 2128.968
## p = 6 q = 1 AIC = 2083.087 BIC = 2113.521
## p = 6 q = 2 AIC = 2084.993 BIC = 2118.471
## p = 6 q = 3 AIC = 2081.777 BIC = 2118.298
## p = 6 q = 4 AIC = 2083.115 BIC = 2122.68
## p = 6 q = 5 AIC = 2084.976 BIC = 2127.584
## p = 6 q = 6 AIC = 2086.338 BIC = 2131.989#ARDL model with p=6 and q=3 gives the best values
gold.model4 <- ardlDlm(x = as.vector(asx$`Gold price`) , y = as.vector(asx$`ASX price`), p = 6 , q = 3 )
summary(gold.model4)##
## Time series regression with "ts" data:
## Start = 7, End = 161
##
## Call:
## dynlm(formula = as.formula(model.text), data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -547.35 -111.67 20.15 114.80 691.34
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 239.21836 95.10550 2.515 0.01299 *
## X.t -1.24010 0.30584 -4.055 8.19e-05 ***
## X.1 1.31981 0.45338 2.911 0.00418 **
## X.2 0.10646 0.46750 0.228 0.82019
## X.3 -0.35071 0.46841 -0.749 0.45525
## X.4 -0.26717 0.45899 -0.582 0.56142
## X.5 0.78380 0.45152 1.736 0.08472 .
## X.6 -0.36216 0.29983 -1.208 0.22908
## Y.1 0.98540 0.08311 11.857 < 2e-16 ***
## Y.2 0.15624 0.11634 1.343 0.18140
## Y.3 -0.18365 0.08272 -2.220 0.02798 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 191.7 on 144 degrees of freedom
## Multiple R-squared: 0.9498, Adjusted R-squared: 0.9464
## F-statistic: 272.7 on 10 and 144 DF, p-value: < 2.2e-16gold.VIF.model4 <- vif(gold.model4$model)
gold.VIF.model4## X.t L(X.t, 1) L(X.t, 2) L(X.t, 3) L(X.t, 4) L(X.t, 5) L(X.t, 6) L(y.t, 1)
## 59.16953 131.52240 141.18835 143.43247 139.07695 135.63989 60.28428 20.34486
## L(y.t, 2) L(y.t, 3)
## 40.73197 21.03241#Sorting the models
gold.aic.sort <- AIC(gold.model$model, gold.model2$model, gold.model4$model)
gold.aic.sort#gold.model4 using ardlDLM has the best AIC score, however it does suffer from multicollinearityDLM Fitting for Oil Price Variable
# Fitting a DLM model
for ( i in 1:6){
oil.model1.1 <- dlm(x = as.vector(asx$`Crude Oil (Brent)_USD/bbl`) , y = as.vector(asx$`ASX price`), q = i)
cat("q = ", i, "AIC = ", AIC(oil.model1.1$model), "BIC = ", BIC(oil.model1.1$model),"\n")
}## q = 1 AIC = 2614.698 BIC = 2626.998
## q = 2 AIC = 2596.715 BIC = 2612.059
## q = 3 AIC = 2579.101 BIC = 2597.477
## q = 4 AIC = 2561.888 BIC = 2583.281
## q = 5 AIC = 2544.936 BIC = 2569.335
## q = 6 AIC = 2527.701 BIC = 2555.091#Best DLM model with lag = 6
oil.model <- dlm( x = as.vector(asx$`Crude Oil (Brent)_USD/bbl`) , y = as.vector(asx$`ASX price`), q = 6)
summary(oil.model)##
## Call:
## lm(formula = model.formula, data = design)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1489.01 -649.84 20.52 631.35 1742.23
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4465.209 189.119 23.611 <2e-16 ***
## x.t 11.039 11.476 0.962 0.338
## x.1 1.035 19.234 0.054 0.957
## x.2 -3.096 19.421 -0.159 0.874
## x.3 4.826 19.469 0.248 0.805
## x.4 -5.191 19.433 -0.267 0.790
## x.5 1.270 19.231 0.066 0.947
## x.6 -4.356 11.386 -0.383 0.703
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 815.1 on 147 degrees of freedom
## Multiple R-squared: 0.07398, Adjusted R-squared: 0.02988
## F-statistic: 1.678 on 7 and 147 DF, p-value: 0.1187
##
## AIC and BIC values for the model:
## AIC BIC
## 1 2527.701 2555.091oil.VIF.model <- vif(oil.model$model)
oil.VIF.model## x.t x.1 x.2 x.3 x.4 x.5 x.6
## 26.21482 74.58090 77.14482 78.64565 79.26209 78.39487 27.75993#Fitting a ploynomial model
for (i in 1:6){
oil.model2.1 <- polyDlm(x = as.vector(asx$`Crude Oil (Brent)_USD/bbl`) , y = as.vector(asx$`ASX price`), q = i , k = i, show.beta = FALSE)
cat("q = ", i, "k = ", i, "AIC = ", AIC(oil.model2.1$model), "BIC = ", BIC(oil.model2.1$model),"\n")
}## q = 1 k = 1 AIC = 2614.698 BIC = 2626.998
## q = 2 k = 2 AIC = 2596.715 BIC = 2612.059
## q = 3 k = 3 AIC = 2579.101 BIC = 2597.477
## q = 4 k = 4 AIC = 2561.888 BIC = 2583.281
## q = 5 k = 5 AIC = 2544.936 BIC = 2569.335
## q = 6 k = 6 AIC = 2527.701 BIC = 2555.091#Best polyDLM model with q = 6 and i = 6
oil.model2 <- polyDlm(x = as.vector(asx$`Crude Oil (Brent)_USD/bbl`) , y = as.vector(asx$`ASX price`), q = 6, k = 6, show.beta = FALSE)
summary(oil.model2)##
## Call:
## "Y ~ (Intercept) + X.t"
##
## Residuals:
## Min 1Q Median 3Q Max
## -1489.01 -649.84 20.52 631.35 1742.23
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4465.2087 189.1191 23.611 <2e-16 ***
## z.t0 11.0386 11.4756 0.962 0.338
## z.t1 56.2691 348.7645 0.161 0.872
## z.t2 -145.3875 662.7956 -0.219 0.827
## z.t3 110.6438 455.2988 0.243 0.808
## z.t4 -36.8002 143.5804 -0.256 0.798
## z.t5 5.5873 21.1137 0.265 0.792
## z.t6 -0.3166 1.1719 -0.270 0.787
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 815.1 on 147 degrees of freedom
## Multiple R-squared: 0.07398, Adjusted R-squared: 0.02988
## F-statistic: 1.678 on 7 and 147 DF, p-value: 0.1187oil.VIF.model2 <- vif(oil.model2$model)
oil.VIF.model2## z.t0 z.t1 z.t2 z.t3 z.t4 z.t5
## 1.220717e+03 1.064677e+07 7.374371e+08 8.288557e+09 2.214560e+10 1.386655e+10
## z.t6
## 1.301240e+09#Applying koyck transformation
oil.model3 <- koyckDlm(x = as.vector(asx$`Crude Oil (Brent)_USD/bbl`) , y = as.vector(asx$`ASX price`))
summary(oil.model3)##
## Call:
## "Y ~ (Intercept) + Y.1 + X.t"
##
## Residuals:
## Min 1Q Median 3Q Max
## -683.91 -108.66 13.68 139.77 762.55
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 209.89536 87.89368 2.388 0.0181 *
## Y.1 0.97537 0.01905 51.193 <2e-16 ***
## X.t -0.99907 0.58045 -1.721 0.0872 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 201.1 on 157 degrees of freedom
## Multiple R-Squared: 0.949, Adjusted R-squared: 0.9483
## Wald test: 1461 on 2 and 157 DF, p-value: < 2.2e-16
##
## alpha beta phi
## Geometric coefficients: 8522.034 -0.9990694 0.9753703oil.VIF.model3 <- vif(oil.model3$model)
oil.VIF.model3## Y.1 X.t
## 1.14038 1.14038#Fitting an ARDL model
for (i in 1:6){
for(j in 1:6){
oil.model4.1 = ardlDlm(x = as.vector(asx$`Crude Oil (Brent)_USD/bbl`) , y = as.vector(asx$`ASX price`), p = i , q = j )
cat("p = ", i, "q = ", j, "AIC = ", AIC(oil.model4.1$model), "BIC = ", BIC(oil.model4.1$model),"\n")
}
}## p = 1 q = 1 AIC = 2146.524 BIC = 2161.9
## p = 1 q = 2 AIC = 2134.107 BIC = 2152.521
## p = 1 q = 3 AIC = 2121.07 BIC = 2142.508
## p = 1 q = 4 AIC = 2109.4 BIC = 2133.85
## p = 1 q = 5 AIC = 2098.335 BIC = 2125.784
## p = 1 q = 6 AIC = 2087.328 BIC = 2117.762
## p = 2 q = 1 AIC = 2132.312 BIC = 2150.726
## p = 2 q = 2 AIC = 2134.235 BIC = 2155.718
## p = 2 q = 3 AIC = 2122.356 BIC = 2146.857
## p = 2 q = 4 AIC = 2110.793 BIC = 2138.299
## p = 2 q = 5 AIC = 2099.752 BIC = 2130.251
## p = 2 q = 6 AIC = 2088.769 BIC = 2122.247
## p = 3 q = 1 AIC = 2121.919 BIC = 2143.357
## p = 3 q = 2 AIC = 2123.835 BIC = 2148.335
## p = 3 q = 3 AIC = 2124.324 BIC = 2151.887
## p = 3 q = 4 AIC = 2112.401 BIC = 2142.963
## p = 3 q = 5 AIC = 2101.35 BIC = 2134.899
## p = 3 q = 6 AIC = 2090.319 BIC = 2126.84
## p = 4 q = 1 AIC = 2111.383 BIC = 2135.832
## p = 4 q = 2 AIC = 2113.294 BIC = 2140.8
## p = 4 q = 3 AIC = 2113.805 BIC = 2144.367
## p = 4 q = 4 AIC = 2114.384 BIC = 2148.003
## p = 4 q = 5 AIC = 2103.342 BIC = 2139.94
## p = 4 q = 6 AIC = 2092.313 BIC = 2131.878
## p = 5 q = 1 AIC = 2097.076 BIC = 2124.525
## p = 5 q = 2 AIC = 2099.041 BIC = 2129.54
## p = 5 q = 3 AIC = 2099.518 BIC = 2133.066
## p = 5 q = 4 AIC = 2099.845 BIC = 2136.443
## p = 5 q = 5 AIC = 2100.917 BIC = 2140.566
## p = 5 q = 6 AIC = 2088.774 BIC = 2131.382
## p = 6 q = 1 AIC = 2086.125 BIC = 2116.559
## p = 6 q = 2 AIC = 2088.117 BIC = 2121.595
## p = 6 q = 3 AIC = 2088.596 BIC = 2125.117
## p = 6 q = 4 AIC = 2088.856 BIC = 2128.421
## p = 6 q = 5 AIC = 2090.082 BIC = 2132.69
## p = 6 q = 6 AIC = 2090.719 BIC = 2136.37#ARDL model with p=6 and q=1 gives the best values
oil.model4 <- ardlDlm(x = as.vector(asx$`Crude Oil (Brent)_USD/bbl`) , y = as.vector(asx$`ASX price`), p = 6 , q = 1 )
summary(oil.model4)##
## Time series regression with "ts" data:
## Start = 7, End = 161
##
## Call:
## dynlm(formula = as.formula(model.text), data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -640.89 -119.93 -6.81 128.21 705.30
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 253.70556 97.08067 2.613 0.0099 **
## X.t 5.62158 2.75536 2.040 0.0431 *
## X.1 -3.06010 4.61525 -0.663 0.5083
## X.2 -3.56288 4.65944 -0.765 0.4457
## X.3 2.93548 4.67105 0.628 0.5307
## X.4 -8.35374 4.66257 -1.792 0.0753 .
## X.5 7.08966 4.61542 1.536 0.1267
## X.6 -1.53920 2.73216 -0.563 0.5741
## Y.1 0.96471 0.01966 49.071 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 195.6 on 146 degrees of freedom
## Multiple R-squared: 0.9471, Adjusted R-squared: 0.9442
## F-statistic: 326.5 on 8 and 146 DF, p-value: < 2.2e-16oil.VIF.model4 <- vif(oil.model4$model)
oil.VIF.model4## X.t L(X.t, 1) L(X.t, 2) L(X.t, 3) L(X.t, 4) L(X.t, 5) L(X.t, 6) L(y.t, 1)
## 26.256967 74.605285 77.145142 78.651000 79.277240 78.446668 27.772195 1.093709#Sorting the models
oil.aic.sort <- AIC(oil.model$model, oil.model2$model, oil.model4$model)
oil.aic.sort#oil.model4 using ardlDLM has the best AIC score, however it does suffer from multicollinearityDLM Fitting for Copper Price Variable
# Fitting a DLM model
for ( i in 1:6){
copper.model1.1 <- dlm(x = as.vector(asx$`Copper_USD/tonne`) , y = as.vector(asx$`ASX price`), q = i)
cat("q = ", i, "AIC = ", AIC(copper.model1.1$model), "BIC = ", BIC(copper.model1.1$model),"\n")
}## q = 1 AIC = 2574.488 BIC = 2586.789
## q = 2 AIC = 2559.356 BIC = 2574.7
## q = 3 AIC = 2544.155 BIC = 2562.531
## q = 4 AIC = 2528.895 BIC = 2550.289
## q = 5 AIC = 2513.265 BIC = 2537.664
## q = 6 AIC = 2497.775 BIC = 2525.166#Best DLM model with lag = 6
copper.model <- dlm( x = as.vector(asx$`Copper_USD/tonne`) , y = as.vector(asx$`ASX price`), q = 6)
summary(copper.model)##
## Call:
## lm(formula = model.formula, data = design)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1200.0 -691.5 -106.3 631.6 1564.7
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.722e+03 1.980e+02 18.799 <2e-16 ***
## x.t 1.582e-01 1.325e-01 1.195 0.234
## x.1 4.879e-02 2.127e-01 0.229 0.819
## x.2 1.864e-02 2.128e-01 0.088 0.930
## x.3 3.026e-02 2.129e-01 0.142 0.887
## x.4 2.216e-02 2.130e-01 0.104 0.917
## x.5 2.903e-03 2.127e-01 0.014 0.989
## x.6 -9.299e-02 1.304e-01 -0.713 0.477
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 740.1 on 147 degrees of freedom
## Multiple R-squared: 0.2366, Adjusted R-squared: 0.2002
## F-statistic: 6.507 on 7 and 147 DF, p-value: 1.086e-06
##
## AIC and BIC values for the model:
## AIC BIC
## 1 2497.775 2525.166copper.VIF.model <- vif(copper.model$model)
copper.VIF.model## x.t x.1 x.2 x.3 x.4 x.5 x.6
## 19.52759 51.96694 53.67132 55.39695 57.01562 58.45737 22.56019#Fitting a ploynomial model
for (i in 1:6){
copper.model2.1 <- polyDlm(x = as.vector(asx$`Copper_USD/tonne`) , y = as.vector(asx$`ASX price`), q = i , k = i, show.beta = FALSE)
cat("q = ", i, "k = ", i, "AIC = ", AIC(copper.model2.1$model), "BIC = ", BIC(copper.model2.1$model),"\n")
}## q = 1 k = 1 AIC = 2574.488 BIC = 2586.789
## q = 2 k = 2 AIC = 2559.356 BIC = 2574.7
## q = 3 k = 3 AIC = 2544.155 BIC = 2562.531
## q = 4 k = 4 AIC = 2528.895 BIC = 2550.289
## q = 5 k = 5 AIC = 2513.265 BIC = 2537.664
## q = 6 k = 6 AIC = 2497.775 BIC = 2525.166#Best polyDLM model with q = 6 and i = 6
copper.model2 <- polyDlm(x = as.vector(asx$`Copper_USD/tonne`) , y = as.vector(asx$`ASX price`), q = 6, k = 6, show.beta = FALSE)
summary(copper.model2)##
## Call:
## "Y ~ (Intercept) + X.t"
##
## Residuals:
## Min 1Q Median 3Q Max
## -1200.0 -691.5 -106.3 631.6 1564.7
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.722e+03 1.980e+02 18.799 <2e-16 ***
## z.t0 1.582e-01 1.325e-01 1.195 0.234
## z.t1 -9.714e-02 3.859e+00 -0.025 0.980
## z.t2 -8.233e-02 7.320e+00 -0.011 0.991
## z.t3 1.015e-01 5.025e+00 0.020 0.984
## z.t4 -3.694e-02 1.584e+00 -0.023 0.981
## z.t5 5.743e-03 2.329e-01 0.025 0.980
## z.t6 -3.307e-04 1.293e-02 -0.026 0.980
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 740.1 on 147 degrees of freedom
## Multiple R-squared: 0.2366, Adjusted R-squared: 0.2002
## F-statistic: 6.507 on 7 and 147 DF, p-value: 1.086e-06copper.VIF.model2 <- vif(copper.model2$model)
copper.VIF.model2## z.t0 z.t1 z.t2 z.t3 z.t4 z.t5
## 9.556226e+02 7.834859e+06 5.452921e+08 6.152844e+09 1.648793e+10 1.034790e+10
## z.t6
## 9.728588e+08#Applying koyck transformation
copper.model3 <- koyckDlm(x = as.vector(asx$`Copper_USD/tonne`) , y = as.vector(asx$`ASX price`))
summary(copper.model3)##
## Call:
## "Y ~ (Intercept) + Y.1 + X.t"
##
## Residuals:
## Min 1Q Median 3Q Max
## -689.64 -108.62 12.78 140.20 771.79
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 189.368812 87.644648 2.161 0.0322 *
## Y.1 0.971621 0.021895 44.376 <2e-16 ***
## X.t -0.005864 0.009517 -0.616 0.5387
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 201.9 on 157 degrees of freedom
## Multiple R-Squared: 0.9485, Adjusted R-squared: 0.9479
## Wald test: 1448 on 2 and 157 DF, p-value: < 2.2e-16
##
## alpha beta phi
## Geometric coefficients: 6672.885 -0.005863623 0.9716211copper.VIF.model3 <- vif(copper.model3$model)
copper.VIF.model3## Y.1 X.t
## 1.493966 1.493966#Fitting an ARDL model
for (i in 1:6){
for(j in 1:6){
copper.model4.1 = ardlDlm(x = as.vector(asx$`Copper_USD/tonne`) , y = as.vector(asx$`ASX price`), p = i , q = j )
cat("p = ", i, "q = ", j, "AIC = ", AIC(copper.model4.1$model), "BIC = ", BIC(copper.model4.1$model),"\n")
}
}## p = 1 q = 1 AIC = 2147.741 BIC = 2163.116
## p = 1 q = 2 AIC = 2135.4 BIC = 2153.813
## p = 1 q = 3 AIC = 2121.12 BIC = 2142.558
## p = 1 q = 4 AIC = 2109.759 BIC = 2134.209
## p = 1 q = 5 AIC = 2099.056 BIC = 2126.505
## p = 1 q = 6 AIC = 2088.466 BIC = 2118.9
## p = 2 q = 1 AIC = 2130.043 BIC = 2148.456
## p = 2 q = 2 AIC = 2132.038 BIC = 2153.52
## p = 2 q = 3 AIC = 2119.241 BIC = 2143.741
## p = 2 q = 4 AIC = 2107.649 BIC = 2135.155
## p = 2 q = 5 AIC = 2097.021 BIC = 2127.52
## p = 2 q = 6 AIC = 2086.386 BIC = 2119.864
## p = 3 q = 1 AIC = 2117.307 BIC = 2138.745
## p = 3 q = 2 AIC = 2119.247 BIC = 2143.748
## p = 3 q = 3 AIC = 2119.696 BIC = 2147.259
## p = 3 q = 4 AIC = 2108.537 BIC = 2139.1
## p = 3 q = 5 AIC = 2097.832 BIC = 2131.38
## p = 3 q = 6 AIC = 2087.217 BIC = 2123.738
## p = 4 q = 1 AIC = 2105.916 BIC = 2130.366
## p = 4 q = 2 AIC = 2107.774 BIC = 2135.28
## p = 4 q = 3 AIC = 2108.608 BIC = 2139.17
## p = 4 q = 4 AIC = 2110.085 BIC = 2143.704
## p = 4 q = 5 AIC = 2099.454 BIC = 2136.052
## p = 4 q = 6 AIC = 2088.844 BIC = 2128.408
## p = 5 q = 1 AIC = 2095.118 BIC = 2122.566
## p = 5 q = 2 AIC = 2096.96 BIC = 2127.459
## p = 5 q = 3 AIC = 2097.887 BIC = 2131.436
## p = 5 q = 4 AIC = 2099.497 BIC = 2136.095
## p = 5 q = 5 AIC = 2101.419 BIC = 2141.067
## p = 5 q = 6 AIC = 2090.789 BIC = 2133.397
## p = 6 q = 1 AIC = 2084.49 BIC = 2114.924
## p = 6 q = 2 AIC = 2086.331 BIC = 2119.809
## p = 6 q = 3 AIC = 2087.163 BIC = 2123.684
## p = 6 q = 4 AIC = 2088.704 BIC = 2128.268
## p = 6 q = 5 AIC = 2090.603 BIC = 2133.211
## p = 6 q = 6 AIC = 2092.599 BIC = 2138.25#ARDL model with p=6 and q=1 gives the best values
copper.model4 <- ardlDlm(x = as.vector(asx$`Copper_USD/tonne`) , y = as.vector(asx$`ASX price`), p = 6 , q = 1 )
summary(copper.model4)##
## Time series regression with "ts" data:
## Start = 7, End = 161
##
## Call:
## dynlm(formula = as.formula(model.text), data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -668.21 -105.68 -6.16 132.96 746.47
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 229.599622 94.131573 2.439 0.0159 *
## X.t 0.082186 0.034857 2.358 0.0197 *
## X.1 -0.025479 0.055929 -0.456 0.6494
## X.2 -0.023929 0.055950 -0.428 0.6695
## X.3 -0.009499 0.055976 -0.170 0.8655
## X.4 -0.011879 0.055979 -0.212 0.8322
## X.5 -0.026840 0.055920 -0.480 0.6320
## X.6 0.008350 0.034357 0.243 0.8083
## Y.1 0.964213 0.021659 44.518 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 194.5 on 146 degrees of freedom
## Multiple R-squared: 0.9476, Adjusted R-squared: 0.9447
## F-statistic: 330.1 on 8 and 146 DF, p-value: < 2.2e-16copper.VIF.model4 <- vif(copper.model4$model)
copper.VIF.model4## X.t L(X.t, 1) L(X.t, 2) L(X.t, 3) L(X.t, 4) L(X.t, 5) L(X.t, 6) L(y.t, 1)
## 19.574613 52.013225 53.687005 55.411058 57.026263 58.465718 22.659672 1.341549#Sorting the models
copper.aic.sort <- AIC(copper.model$model, copper.model2$model, copper.model4$model)
copper.aic.sort#copper.model4 using ardlDLM has the best AIC score, however it does suffer from multicollinearity#Comparing the models from all 3 variables
aic.sort <- AIC(gold.model4$model, oil.model4$model, copper.model4$model)
aic.sortgold.VIF.model4## X.t L(X.t, 1) L(X.t, 2) L(X.t, 3) L(X.t, 4) L(X.t, 5) L(X.t, 6) L(y.t, 1)
## 59.16953 131.52240 141.18835 143.43247 139.07695 135.63989 60.28428 20.34486
## L(y.t, 2) L(y.t, 3)
## 40.73197 21.03241Conclusion for best DLM for ASX Price Index
The gold variable modeled using the ardlDLM provides the best DLM for the ASX Price Index due to it having the best R-squared score. Gold also has the lowest AIC score. However, all 3 variable models, including gold, suffer from multicollinearity due to high VIF scores.