Reproduction of House Price Modeling
Abdul Samad
June 26, 2021
Service as a Freelancer
If you want hire me as a freelancer please use the following links.
- For Fiverr https://www.fiverr.com/dataanalyst_85
- For Upwork https://www.upwork.com/freelancers/~017667e5ef4347dcb4
2years Sample Analysis of Monthly Data
# specify the type of data to be used
datasample <- read_excel("hpimonthly.xls")
# specify size of out of sample data
outsampsize <- 60
# specify frequency of timeseries data
fq = 12
# Corresponding graphs will also be generated
datset <- datasample
hpx <- datset # load house price data series
hpx$year <- as.numeric(format(hpx$Date,"%Y"))
hpx$month <- as.numeric(format(hpx$Date, "%m"))
# set the size of out-sample data
outsampsize = outsampsize
#### ARIMA - GARCH model for house price index
###############################################
# prepare full dataset for analysis
N <- nrow(hpx) # length of dataset
M <- N - outsampsize # index for training sample
hpx.main <- hpx$Price[1:N] # main sample for analysis
mainsample <- ts(hpx.main, start = c(hpx$year[1], hpx$month[1]), end = c(hpx$year[N], hpx$month[N]), frequency = fq)
# generate log returns for main sample and holdout
rets_mainsample <- diff(log(mainsample)) # main sample log-returns
rets_holdoutsample <- rets_mainsample[M:(N-1)] # holdout sample log-returns
########################################
# ARMA-EGARCH model
########################################
# specify the model
model.garch = ugarchspec(mean.model = list(armaOrder = c(4,3), include.mean = TRUE),variance.model = list(garchOrder = c(1,1), model = "eGARCH"),
distribution.model = "norm")
# fit the specified model
model.garch.fit = ugarchfit(data = rets_mainsample, spec = model.garch,
out.sample = outsampsize, solver = "hybrid" )
# check coefficients of fitted model
model.garch.fit@fit$robust.matcoef## Estimate Std. Error t value Pr(>|t|)
## mu 0.005871745 6.189227e-06 948.704044 0.00000000
## ar1 0.942194430 4.970087e-04 1895.730220 0.00000000
## ar2 -0.762975366 5.600831e-04 -1362.253831 0.00000000
## ar3 -0.033770149 3.664604e-05 -921.522387 0.00000000
## ar4 0.388518050 2.218478e-04 1751.282263 0.00000000
## ma1 -0.685798639 2.872588e-04 -2387.389273 0.00000000
## ma2 1.033908138 4.850121e-05 21317.161344 0.00000000
## ma3 0.032425864 4.125027e-05 786.076306 0.00000000
## omega -1.134627627 4.953247e-01 -2.290674 0.02198225
## alpha1 -0.126142856 5.750292e-02 -2.193677 0.02825862
## beta1 0.878108776 5.424889e-02 16.186669 0.00000000
## gamma1 0.163384889 1.332465e-01 1.226185 0.22012897library(ggplot2)
library(broom)
library(zoo)
fitted_values <- fitted(model.garch.fit)
resid_values <- residuals(model.garch.fit)
actual_values <- fitted_values[, 1]+resid_values[, 1]
fitted_values <- merge(fitted_values,actual_values)
names(fitted_values) <- c("fitted", "actual")
index(fitted_values) <- as.Date(index(fitted_values))
tidy(fitted_values) %>% ggplot(aes(x = index, y = value, color = series)) + geom_line()
plot(model.garch.fit, which = "all")##
## please wait...calculating quantiles...
# out of sample forecast
library(reshape2)
library(ggplot2)
modelfor = ugarchforecast(model.garch.fit, data = NULL, n.ahead = outsampsize,
n.roll= 0, out.sample = outsampsize)
outsampleforecast <- data.frame(date = hpx$Date[M:(N - 1)],
seriesforcast = modelfor@forecast$seriesFor,
actual = rets_holdoutsample)
rownames(outsampleforecast) <- NULL
colnames(outsampleforecast) <- c('date', 'seriesforcast', 'holdoutsample')
meltedoutsampleforecast <- melt(outsampleforecast, id.vars = 'date')
head(outsampleforecast)## date seriesforcast holdoutsample
## 1 2013-09-01 0.003171730 0.008971460
## 2 2013-10-01 0.007305870 0.005102953
## 3 2013-11-01 0.006319372 0.007186950
## 4 2013-12-01 0.006666740 0.003777293
## 5 2014-01-01 0.005181821 0.007647984
## 6 2014-02-01 0.005157208 0.013506314###############################################################################
# Fitting GBM
###############################################################################
X <- ts(hpx.main, start = c(hpx$year[1], hpx$month[1]), end = c(hpx$year[M], hpx$month[M]), frequency = fq)
#' redefine main sample as X
dcBS <- function(x, t, x0 , theta , log = TRUE ){
ml <- log(x0) + (theta[1] - theta[2]^2/2)*t
sl <- sqrt(t) * theta[2]
lik <- dlnorm(x, meanlog = ml , sdlog = sl , log = TRUE )
if(!log)
lik <- exp(lik)
lik
}
BS.lik <- function(theta1 , theta2){
n <- length(X)
dt <- deltat(X)
-sum(dcBS(x = X[2:n], t = dt, x0 = X[1:(n - 1)], theta = c(theta1, theta2), log = TRUE ))
}
gbm_mle <- mle(BS.lik, start=list(theta1=.01, theta2=.01),
method="L-BFGS-B",lower=c(0.001,0.001), upper=c(1,1))
# maximum likelihood estimates
coef(gbm_mle)## theta1 theta2
## 0.05277341 0.04113015# Method of Moments Estimation
# quick test with plug-in estimators
returns_X <- diff(log(X)) # main sample log-returns
alpha <- mean(returns_X)/(1/fq)
sigma <- sqrt(var(returns_X)/(1/fq))
mu <- alpha + 0.5*sigma^2
gbm_mm <- cbind(theta1 = mu, theta2 = sigma)
gbm_mm## theta1 theta2
## [1,] 0.05277335 0.0411939# GMM Estimation
DELTA = 1/fq
u <- function(x, y, theta, DELTA){
c.mean <- theta[1]*DELTA
c.var <- theta[2]^2*DELTA
cbind(x - c.mean,y*(x - c.mean), c.var - (x - c.mean)^2,
y*(c.var - (x - c.mean)^2))
} # gmm function definition
gbm_gmm <- sde::gmm(returns_X, u, guess = as.numeric(coef(gbm_mle)), lower = c(0.01, 0.01), upper = c(1,1), tol1 = 1e-3, tol2 = 1e-3, maxiter = 1000 );##
## Dimension of parameter space set to 2
## Initial values for the optimization algorithm
## [1] 0.05277341 0.04113015
##
## Optimization contraints
## lower upper
## theta1 0.01 1
## theta2 0.01 1
##
## Running optimizer...
##
## First stage estimates:
## theta1 theta2
## 0.05275919 0.04113015
##
## Starting second stage...
## theta1 theta2 Q1 |theta1-theta2|
## 0.046657025 0.041200557 0.004973025 0.006172567
## theta1 theta2 Q1 |theta1-theta2|
## 0.041429050 0.041238387 0.004794165 0.005265805
## theta1 theta2 Q1 |theta1-theta2|
## 0.037125185 0.041259777 0.004672235 0.004325254
## theta1 theta2 Q1 |theta1-theta2|
## 0.033800626 0.041238805 0.004611018 0.003345531
## theta1 theta2 Q1 |theta1-theta2|
## 0.031156062 0.041209222 0.004591138 0.002674147
## theta1 theta2 Q1 |theta1-theta2|
## 0.029112696 0.041164160 0.004594039 0.002088428
## theta1 theta2 Q1 |theta1-theta2|
## 0.027519783 0.041096136 0.004607621 0.001660938
## theta1 theta2 Q1 |theta1-theta2|
## 0.026241371 0.041024655 0.004624570 0.001349893
## theta1 theta2 Q1 |theta1-theta2|
## 0.025215932 0.040951358 0.004642234 0.001098736
## theta1 theta2 Q1 |theta1-theta2|
## 0.0243806559 0.0408811404 0.0046588144 0.0009054935gbm_gmm$par## theta1 theta2
## 0.02438066 0.04088114# gather the estimated gbm coefficients from the parameter estimation methods i.e. MLE, GMM & MM
gbm.est <- data.frame(matrix(unlist(list(MLE = coef(gbm_mle), GMM = gbm_gmm$par, MM = gbm_mm )), nrow=3, byrow=T))
rownames(gbm.est) <- c("MLE", "GMM", "MM")
colnames(gbm.est) <- c("mu", "sigma")
#######################################################################################################
#### Forecast under each model
#######################################################################################################
# Forecast the GBM outsampsize-months ahead based on estimation method
sim_mle <- snssde1d(drift = expression(gbm.est[1,1]*x), diffusion = expression(gbm.est[1,2]*x),
M = 1, N = outsampsize, x0 = rets_mainsample[M], Dt = 1/fq)$X[-1]
sim_gmm <- snssde1d(drift = expression(gbm.est[2,1]*x), diffusion = expression(gbm.est[2,2]*x),
M = 1, N = outsampsize, x0 = rets_mainsample[M], Dt = 1/fq)$X[-1]
sim_mm <- snssde1d(drift = expression(gbm.est[3,1]*x), diffusion = expression(gbm.est[3,2]*x),
M = 1, N = outsampsize, x0 = rets_mainsample[M], Dt = 1/fq)$X[-1]
#######################################################################################################
#### Pull all forecasted models together for Dielbod-Mariano test of predictive accuracy
#######################################################################################################
outsampleforecast$gbm_mle <- sim_mle
outsampleforecast$gbm_mm <- sim_mm
outsampleforecast$gbm_gmm <- sim_gmm
forcast_errors <- data.frame(date = outsampleforecast$date,
MLE = (outsampleforecast$holdoutsample - outsampleforecast$gbm_mle),
GMM = (outsampleforecast$holdoutsample - outsampleforecast$gbm_gmm),
MM = (outsampleforecast$holdoutsample - outsampleforecast$gbm_mm),
ARMA_EGARCH = (outsampleforecast$holdoutsample - outsampleforecast$seriesforcast))
head(forcast_errors)## date MLE GMM MM ARMA_EGARCH
## 1 2013-09-01 -9.273274e-05 -9.508906e-05 6.863196e-05 0.005799730
## 2 2013-10-01 -3.883775e-03 -3.951178e-03 -3.860936e-03 -0.002202917
## 3 2013-11-01 -1.763799e-03 -1.761759e-03 -1.932592e-03 0.000867578
## 4 2013-12-01 -5.085726e-03 -5.349017e-03 -5.297739e-03 -0.002889447
## 5 2014-01-01 -1.162558e-03 -1.414556e-03 -1.415947e-03 0.002466164
## 6 2014-02-01 4.770937e-03 4.472365e-03 4.270152e-03 0.008349106meltedforecasterrors <- melt(forcast_errors, id.vars = 'date')
ggplot(data = meltedforecasterrors, aes(x = date, y = value, color = variable)) + geom_line()
head(outsampleforecast)## date seriesforcast holdoutsample gbm_mle gbm_mm gbm_gmm
## 1 2013-09-01 0.003171730 0.008971460 0.009064193 0.008902828 0.009066549
## 2 2013-10-01 0.007305870 0.005102953 0.008986729 0.008963889 0.009054131
## 3 2013-11-01 0.006319372 0.007186950 0.008950749 0.009119542 0.008948709
## 4 2013-12-01 0.006666740 0.003777293 0.008863020 0.009075033 0.009126311
## 5 2014-01-01 0.005181821 0.007647984 0.008810543 0.009063932 0.009062541
## 6 2014-02-01 0.005157208 0.013506314 0.008735376 0.009236162 0.009033949# Estimate Root Mean Squared Error (RMSE)
library(Metrics)
tab_rmse <- data.frame(rbind(GBM.MLE = rmse(outsampleforecast$holdoutsample,outsampleforecast$gbm_mle),
GBM.GMM = rmse(outsampleforecast$holdoutsample,outsampleforecast$gbm_gmm),
GBM.MM = rmse(outsampleforecast$holdoutsample,outsampleforecast$gbm_mm),
ARMA_EGARCH = rmse(outsampleforecast$holdoutsample,outsampleforecast$seriesforcast)))
colnames(tab_rmse) <- c('RMSE')
tab_rmse## RMSE
## GBM.MLE 0.007977198
## GBM.GMM 0.008480500
## GBM.MM 0.009083884
## ARMA_EGARCH 0.006528992# Estimate Mean Absolute Error (MAE)
tab_mae <- data.frame(rbind(GBM.MLE = mae(outsampleforecast$holdoutsample,outsampleforecast$gbm_mle),
GBM.GMM = mae(outsampleforecast$holdoutsample,outsampleforecast$gbm_gmm),
GBM.MM = mae(outsampleforecast$holdoutsample,outsampleforecast$gbm_mm),
ARMA_GARCH = mae(outsampleforecast$holdoutsample,outsampleforecast$seriesforcast)))
colnames(tab_mae) <- c('MAE')
tab_mae## MAE
## GBM.MLE 0.006794519
## GBM.GMM 0.007289831
## GBM.MM 0.007855656
## ARMA_GARCH 0.005408543rmse_mae <- cbind(tab_rmse, tab_mae)
rmse_mae## RMSE MAE
## GBM.MLE 0.007977198 0.006794519
## GBM.GMM 0.008480500 0.007289831
## GBM.MM 0.009083884 0.007855656
## ARMA_EGARCH 0.006528992 0.005408543# Diebold & Mariano Test (dmt)
library(forecast)
# estimate test statistic of the dm_test
dm_stats <- data.frame(rbind(MLE_GMM.statistic = dm.test(forcast_errors$MLE,forcast_errors$GMM, alternative = "two.sided", h = 1, power = 2)$statistic,
MLE_MM.statistic = dm.test(forcast_errors$MLE,forcast_errors$MM, alternative = "two.sided", h = 1, power = 2)$statistic,
MLE_ARMAEGARCH.statistic = dm.test(forcast_errors$MLE,forcast_errors$ARMA_EGARCH, alternative = "two.sided", h = 1, power = 2)$statistic,
GMM_MM.statistic = dm.test(forcast_errors$GMM,forcast_errors$MM, alternative = "two.sided", h = 1, power = 2)$statistic,
GMM_ARMAEGARCH.statistic = dm.test(forcast_errors$GMM,forcast_errors$ARMA_EGARCH, alternative = "two.sided", h = 1, power = 2)$statistic,
MM_ARMAEGARCH.statistic = dm.test(forcast_errors$MM,forcast_errors$ARMA_EGARCH, alternative = "two.sided", h = 1, power = 2)$statistic))
colnames(dm_stats) <- c('Statistic')
# estimate corresponding p_value
dm_pvalue <- data.frame(rbind(MLE_GMM.pvalue = dm.test(forcast_errors$MLE,forcast_errors$GMM, alternative = "two.sided", h = 1, power = 2)$p.value,
MLE_MM.pvalue = dm.test(forcast_errors$MLE,forcast_errors$MM, alternative = "two.sided", h = 1, power = 2)$p.value,
MLE_ARMAEGARCH.pvalue = dm.test(forcast_errors$MLE,forcast_errors$ARMA_EGARCH, alternative = "two.sided", h = 1, power = 2)$p.value,
GMM_MM.pvalue = dm.test(forcast_errors$GMM,forcast_errors$MM, alternative = "two.sided", h = 1, power = 2)$p.value,
GMM_ARMAEGARCH.pvalue = dm.test(forcast_errors$GMM,forcast_errors$ARMA_EGARCH, alternative = "two.sided", h = 1, power = 2)$p.value,
MM_ARMAEGARCH.pvalue = dm.test(forcast_errors$MM,forcast_errors$ARMA_EGARCH, alternative = "two.sided", h = 1, power = 2)$p.value))
colnames(dm_pvalue) <- c('P-Value')
dmtest <- cbind(dm_stats, dm_pvalue)
dmtest## Statistic P-Value
## MLE_GMM.statistic -5.700811 4.033272e-07
## MLE_MM.statistic -6.735876 7.562238e-09
## MLE_ARMAEGARCH.statistic 4.520558 3.024054e-05
## GMM_MM.statistic -6.597257 1.295233e-08
## GMM_ARMAEGARCH.statistic 4.992343 5.602392e-06
## MM_ARMAEGARCH.statistic 5.537808 7.452383e-07results <- list(ArmaEgarch.parameters = model.garch.fit@fit$robust.matcoef, Gbm.parameters = gbm.est,
MAE.RMSE = rmse_mae, Diebold.Mariano = dmtest)
# output of analysis
results## $ArmaEgarch.parameters
## Estimate Std. Error t value Pr(>|t|)
## mu 0.005871745 6.189227e-06 948.704044 0.00000000
## ar1 0.942194430 4.970087e-04 1895.730220 0.00000000
## ar2 -0.762975366 5.600831e-04 -1362.253831 0.00000000
## ar3 -0.033770149 3.664604e-05 -921.522387 0.00000000
## ar4 0.388518050 2.218478e-04 1751.282263 0.00000000
## ma1 -0.685798639 2.872588e-04 -2387.389273 0.00000000
## ma2 1.033908138 4.850121e-05 21317.161344 0.00000000
## ma3 0.032425864 4.125027e-05 786.076306 0.00000000
## omega -1.134627627 4.953247e-01 -2.290674 0.02198225
## alpha1 -0.126142856 5.750292e-02 -2.193677 0.02825862
## beta1 0.878108776 5.424889e-02 16.186669 0.00000000
## gamma1 0.163384889 1.332465e-01 1.226185 0.22012897
##
## $Gbm.parameters
## mu sigma
## MLE 0.05277341 0.04113015
## GMM 0.02438066 0.04088114
## MM 0.05277335 0.04119390
##
## $MAE.RMSE
## RMSE MAE
## GBM.MLE 0.007977198 0.006794519
## GBM.GMM 0.008480500 0.007289831
## GBM.MM 0.009083884 0.007855656
## ARMA_EGARCH 0.006528992 0.005408543
##
## $Diebold.Mariano
## Statistic P-Value
## MLE_GMM.statistic -5.700811 4.033272e-07
## MLE_MM.statistic -6.735876 7.562238e-09
## MLE_ARMAEGARCH.statistic 4.520558 3.024054e-05
## GMM_MM.statistic -6.597257 1.295233e-08
## GMM_ARMAEGARCH.statistic 4.992343 5.602392e-06
## MM_ARMAEGARCH.statistic 5.537808 7.452383e-075years Sample Analysis of Monthly Data
# specify the type of data to be used
datasample <- read_excel("hpiquarterly.xlsx")
# specify size of out of sample data
outsampsize <- 20
# specify frequency of timeseries data
fq = 4 # (NOTE: 12 = monthly, 4 = quarterly)
# Corresponding graphs will also be generated
datset <- datasample
hpx <- datset # load house price data series
hpx$year <- as.numeric(format(hpx$Date,"%Y"))
hpx$month <- as.numeric(format(hpx$Date, "%m"))
# set the size of out-sample data
outsampsize = outsampsize
#### ARIMA - GARCH model for house price index
###############################################
# prepare full dataset for analysis
# prepare full dataset for analysis
N <- nrow(hpx) # length of dataset
M <- N - outsampsize # index for training sample
hpx.main <- hpx$Price[1:N] # main sample for analysis
mainsample <- ts(hpx.main, start = c(hpx$year[1], hpx$month[1]), end = c(hpx$year[N], hpx$month[N]), frequency = fq)
# generate log returns for main sample and holdout
rets_mainsample <- diff(log(mainsample)) # main sample log-returns
rets_holdoutsample <- rets_mainsample[M:(N-1)] # holdout sample log-returns
########################################
# ARMA-EGARCH model
########################################
# specify the model
model.garch = ugarchspec(mean.model = list(armaOrder = c(4,3), include.mean = TRUE),
variance.model = list(garchOrder = c(1,1), model = "eGARCH"),
distribution.model = "norm")
# fit the specified model
model.garch.fit = ugarchfit(data = rets_mainsample, spec = model.garch, out.sample = outsampsize, solver = "hybrid" )## Warning in arima(data, order = c(modelinc[2], 0, modelinc[3]), include.mean =
## modelinc[1], : possible convergence problem: optim gave code = 1# check coefficients of fitted model
model.garch.fit@fit$robust.matcoef## Estimate Std. Error t value Pr(>|t|)
## mu -0.002373322 1.529654e-06 -1551.5413278 0.000000e+00
## ar1 -0.091661654 1.822190e-03 -50.3030259 0.000000e+00
## ar2 -0.119877071 8.506420e-04 -140.9254029 0.000000e+00
## ar3 -0.108157906 2.352414e-03 -45.9774163 0.000000e+00
## ar4 0.891410859 7.415993e-03 120.2011514 0.000000e+00
## ma1 0.990996126 1.961266e-04 5052.8388534 0.000000e+00
## ma2 0.952121683 4.419577e-04 2154.3276762 0.000000e+00
## ma3 1.040280467 3.612554e-04 2879.6260846 0.000000e+00
## omega -3.583975692 1.138476e+00 -3.1480476 1.643649e-03
## alpha1 0.076016456 7.675645e-02 0.9903592 3.219986e-01
## beta1 0.574464515 1.382349e-01 4.1557118 3.242764e-05
## gamma1 0.321089259 1.480378e-01 2.1689681 3.008511e-02#plot(model.garch.fit, which = "all")
fitted_values <- fitted(model.garch.fit)
resid_values <- residuals(model.garch.fit)
actual_values <- fitted_values[, 1]+resid_values[, 1]
fitted_values <- merge(fitted_values,actual_values)
names(fitted_values) <- c("fitted", "actual")
index(fitted_values) <- as.Date(index(fitted_values))
tidy(fitted_values) %>% ggplot(aes(x = index, y = value, color = series)) + geom_line()
plot(model.garch.fit, which = "all")##
## please wait...calculating quantiles...
# out of sample forecast
modelfor = ugarchforecast(model.garch.fit, data = NULL, n.ahead = outsampsize, n.roll
= 0, out.sample = outsampsize)
outsampleforecast <- data.frame(date = hpx$Date[M:(N - 1)],
seriesforcast = modelfor@forecast$seriesFor,
actual = rets_holdoutsample)
rownames(outsampleforecast) <- NULL
colnames(outsampleforecast) <- c('date', 'seriesforcast', 'holdoutsample')
meltedoutsampleforecast <- melt(outsampleforecast, id.vars = 'date')
ggplot(data = meltedoutsampleforecast, aes(x = date, y = value, color = variable)) + geom_line()
head(outsampleforecast)## date seriesforcast holdoutsample
## 1 2013-10-01 0.014950458 0.020873837
## 2 2014-01-01 0.041556938 0.046189182
## 3 2014-04-01 0.021486587 0.012073571
## 4 2014-07-01 0.008619872 0.001012503
## 5 2014-10-01 0.004449968 -0.002309500
## 6 2015-01-01 0.032262685 0.029741852###############################################################################
# Fitting GBM
###############################################################################
fx <- expression(theta[1]*x) # the drift coefficient for house price gbm model
gx <- expression(theta[2]*x) # the diffusion coefficient of the house price gbm model
gbm_mle <- fitsde(data = mainsample, drift = fx, diffusion = gx,
start = list(theta1 = 0.05,theta2 = 0.05), pmle = "euler") # method of estimation
gbm_mle##
## Call:
## fitsde(data = mainsample, drift = fx, diffusion = gx, start = list(theta1 = 0.05,
## theta2 = 0.05), pmle = "euler")
##
## Coefficients:
## theta1 theta2
## 0.07352602 0.05006677# Method of Moments Estimation
# quick test with plug-in estimators
alpha <- mean(rets_mainsample)/(1/fq)
sigma <- sqrt(var(rets_mainsample)/(1/fq))
mu <- alpha + 0.5*sigma^2
gbm_mm <- list(MM = cbind(theta1 = mu, theta2 = sigma))
gbm_mm## $MM
## theta1 theta2
## [1,] 0.07285854 0.04887302# GMM Estimation
DELTA = 1/fq
u <- function(x, y, theta, DELTA){
c.mean <- theta[1] * DELTA
c.var <- theta[2]^2 * DELTA
cbind(x - c.mean,y*(x - c.mean), c.var - (x - c.mean)^2,
y*(c.var - (x - c.mean)^2))
} # gmm function definition
gbm_gmm <- sde::gmm(rets_mainsample, u, guess = as.numeric(coef(gbm_mle)), lower = c(0.01, 0.01), upper = c(1,1), tol1 = 1e-4, tol2 = 1e-4, maxiter = 1000 )##
## Dimension of parameter space set to 2
## Initial values for the optimization algorithm
## [1] 0.07352602 0.05006677
##
## Optimization contraints
## lower upper
## theta1 0.01 1
## theta2 0.01 1
##
## Running optimizer...
##
## First stage estimates:
## theta1 theta2
## 0.07196420 0.04978743
##
## Starting second stage...
## theta1 theta2 Q1 |theta1-theta2|
## 0.068189946 0.048637218 0.005214742 0.004924470
## theta1 theta2 Q1 |theta1-theta2|
## 0.064695062 0.047484827 0.005181181 0.004647276
## theta1 theta2 Q1 |theta1-theta2|
## 0.061495129 0.046334803 0.005152396 0.004349957
## theta1 theta2 Q1 |theta1-theta2|
## 0.058581856 0.045190977 0.005128344 0.004057099
## theta1 theta2 Q1 |theta1-theta2|
## 0.055947928 0.044056677 0.005108665 0.003768227
## theta1 theta2 Q1 |theta1-theta2|
## 0.053582019 0.042934906 0.005092927 0.003487681
## theta1 theta2 Q1 |theta1-theta2|
## 0.051468643 0.041828325 0.005080635 0.003219957
## theta1 theta2 Q1 |theta1-theta2|
## 0.049590212 0.040739367 0.005071276 0.002967389
## theta1 theta2 Q1 |theta1-theta2|
## 0.047930028 0.039670388 0.005064360 0.002729163
## theta1 theta2 Q1 |theta1-theta2|
## 0.046461303 0.038621475 0.005059479 0.002517639
## theta1 theta2 Q1 |theta1-theta2|
## 0.045163097 0.037595704 0.005056157 0.002323977
## theta1 theta2 Q1 |theta1-theta2|
## 0.044027959 0.036593680 0.005054090 0.002137161
## theta1 theta2 Q1 |theta1-theta2|
## 0.043036528 0.035613952 0.005053111 0.001971158
## theta1 theta2 Q1 |theta1-theta2|
## 0.042171860 0.034659119 0.005052863 0.001819501
## theta1 theta2 Q1 |theta1-theta2|
## 0.041420025 0.033729678 0.005053208 0.001681276
## theta1 theta2 Q1 |theta1-theta2|
## 0.040768463 0.032825977 0.005054000 0.001555263
## theta1 theta2 Q1 |theta1-theta2|
## 0.040205921 0.031948231 0.005055126 0.001440288
## theta1 theta2 Q1 |theta1-theta2|
## 0.039722366 0.031096550 0.005056495 0.001335236
## theta1 theta2 Q1 |theta1-theta2|
## 0.039309856 0.030271083 0.005058038 0.001237977
## theta1 theta2 Q1 |theta1-theta2|
## 0.038959597 0.029471627 0.005059711 0.001149715
## theta1 theta2 Q1 |theta1-theta2|
## 0.038663092 0.028697851 0.005061476 0.001070281
## theta1 theta2 Q1 |theta1-theta2|
## 0.0384154931 0.0279496830 0.0050632948 0.0009957665
## theta1 theta2 Q1 |theta1-theta2|
## 0.0382112955 0.0272268655 0.0050651534 0.0009270152
## theta1 theta2 Q1 |theta1-theta2|
## 0.0380456171 0.0265291024 0.0050670379 0.0008634415
## theta1 theta2 Q1 |theta1-theta2|
## 0.0379141142 0.0258560669 0.0050689390 0.0008045384
## theta1 theta2 Q1 |theta1-theta2|
## 0.037812907 0.025207409 0.005070850 0.000749865
## theta1 theta2 Q1 |theta1-theta2|
## 0.0377385214 0.0245827686 0.0050727666 0.0006990261
## theta1 theta2 Q1 |theta1-theta2|
## 0.0376878064 0.0239817785 0.0050746826 0.0006517051
## theta1 theta2 Q1 |theta1-theta2|
## 0.0376578955 0.0234040698 0.0050765928 0.0006076196
## theta1 theta2 Q1 |theta1-theta2|
## 0.0376461577 0.0228492758 0.0050784901 0.0005665318
## theta1 theta2 Q1 |theta1-theta2|
## 0.0376219019 0.0223265734 0.0050784561 0.0005469582
## theta1 theta2 Q1 |theta1-theta2|
## 0.0376422736 0.0218153727 0.0050821754 0.0005315724
## theta1 theta2 Q1 |theta1-theta2|
## 0.0376732337 0.0213266072 0.0050838678 0.0005197256
## theta1 theta2 Q1 |theta1-theta2|
## 0.0377131863 0.0208628249 0.0050849079 0.0005037348
## theta1 theta2 Q1 |theta1-theta2|
## 0.0377627063 0.0204228912 0.0050860923 0.0004894537
## theta1 theta2 Q1 |theta1-theta2|
## 0.0378198338 0.0200037552 0.0050876485 0.0004762635
## theta1 theta2 Q1 |theta1-theta2|
## 0.037881543 0.019603977 0.005089244 0.000461488
## theta1 theta2 Q1 |theta1-theta2|
## 0.0379456185 0.0192232140 0.0050907149 0.0004448379
## theta1 theta2 Q1 |theta1-theta2|
## 0.0378904524 0.0189038064 0.0050830831 0.0003745737
## theta1 theta2 Q1 |theta1-theta2|
## 0.037815148 0.018591932 0.005086611 0.000387179
## theta1 theta2 Q1 |theta1-theta2|
## 0.0377420344 0.0182964366 0.0050872051 0.0003686088
## theta1 theta2 Q1 |theta1-theta2|
## 0.0376795009 0.0180140585 0.0050878092 0.0003449116
## theta1 theta2 Q1 |theta1-theta2|
## 0.0376599806 0.0177555389 0.0050845105 0.0002780399
## theta1 theta2 Q1 |theta1-theta2|
## 0.0376482882 0.0175069649 0.0050852961 0.0002602664
## theta1 theta2 Q1 |theta1-theta2|
## 0.0376434682 0.0172677593 0.0050861716 0.0002440256
## theta1 theta2 Q1 |theta1-theta2|
## 0.0376455858 0.0170376531 0.0050870906 0.0002322239
## theta1 theta2 Q1 |theta1-theta2|
## 0.0375806230 0.0168236039 0.0050911058 0.0002790119
## theta1 theta2 Q1 |theta1-theta2|
## 0.0375587050 0.0166254028 0.0050871661 0.0002201191
## theta1 theta2 Q1 |theta1-theta2|
## 0.037537707 0.016420834 0.005092083 0.000225567
## theta1 theta2 Q1 |theta1-theta2|
## 0.0374926236 0.0162597628 0.0050866397 0.0002061544
## theta1 theta2 Q1 |theta1-theta2|
## 0.0374405709 0.0161127463 0.0050867028 0.0001990692
## theta1 theta2 Q1 |theta1-theta2|
## 0.0373985747 0.0159720237 0.0050868142 0.0001827188
## theta1 theta2 Q1 |theta1-theta2|
## 0.0373603772 0.0158397247 0.0050869187 0.0001704966
## theta1 theta2 Q1 |theta1-theta2|
## 0.0373323593 0.0157125420 0.0050870361 0.0001552006
## theta1 theta2 Q1 |theta1-theta2|
## 0.0373330389 0.0155882638 0.0050869312 0.0001249578
## theta1 theta2 Q1 |theta1-theta2|
## 0.0373203400 0.0154718473 0.0050873365 0.0001291155
## theta1 theta2 Q1 |theta1-theta2|
## 0.0373044897 0.0153650997 0.0050873859 0.0001225978
## theta1 theta2 Q1 |theta1-theta2|
## 0.0372907495 0.0152650015 0.0050874835 0.0001138384
## theta1 theta2 Q1 |theta1-theta2|
## 0.0372706328 0.0151750298 0.0050875532 0.0001100885
## theta1 theta2 Q1 |theta1-theta2|
## 0.0372536972 0.0150901547 0.0050876454 0.0001018108
## theta1 theta2 Q1 |theta1-theta2|
## 3.724100e-02 1.500946e-02 5.087736e-03 9.338707e-05gbm_gmm$par## theta1 theta2
## 0.03724100 0.01500946# gather the estimated gbm coefficients from the parameter estimation methods i.e. MLE, GMM & MM
gbm.est <- data.frame(matrix(unlist(list(MLE = coef(gbm_mle), GMM = gbm_gmm$par, MM = gbm_mm )), nrow=3, byrow=T))
rownames(gbm.est) <- c("MLE", "GMM", "MM")
colnames(gbm.est) <- c("mu", "sigma")
#######################################################################################################
#### Forecast under each model
#######################################################################################################
# Forecast the GBM outsampsize-months ahead based on estimation method
set.seed(111)
sim_mle <- snssde1d(drift = expression(gbm.est[1,1]*x), diffusion = expression(gbm.est[1,2]*x),
M = 1, N = outsampsize, x0 = rets_mainsample[M], Dt = 1/fq)$X[-1]
set.seed(222)
sim_gmm <- snssde1d(drift = expression(gbm.est[2,1]*x), diffusion = expression(gbm.est[2,2]*x),
M = 1, N = outsampsize, x0 = rets_mainsample[M], Dt = 1/fq)$X[-1]
set.seed(333)
sim_mm <- snssde1d(drift = expression(gbm.est[3,1]*x), diffusion = expression(gbm.est[3,2]*x),
M = 1, N = outsampsize, x0 = rets_mainsample[M], Dt = 1/fq)$X[-1]
#######################################################################################################
#### Pull all forecasted models together for Dielbod-Mariano test of predictive accuracy
#######################################################################################################
outsampleforecast$gbm_mle <- sim_mle
outsampleforecast$gbm_mm <- sim_mm
outsampleforecast$gbm_gmm <- sim_gmm
forcast_errors <- data.frame(date = outsampleforecast$date,
MLE = (outsampleforecast$holdoutsample - outsampleforecast$gbm_mle),
GMM = (outsampleforecast$holdoutsample - outsampleforecast$gbm_gmm),
MM = (outsampleforecast$holdoutsample - outsampleforecast$gbm_mm),
ARMA_EGARCH = (outsampleforecast$holdoutsample - outsampleforecast$seriesforcast))
head(forcast_errors)## date MLE GMM MM ARMA_EGARCH
## 1 2013-10-01 -0.0005066054 -0.0004274016 -0.0003379684 0.005923379
## 2 2014-01-01 0.0245927528 0.0246899259 0.0235881839 0.004632244
## 3 2014-04-01 -0.0097513591 -0.0098486708 -0.0098061901 -0.009413016
## 4 2014-07-01 -0.0199557112 -0.0210512873 -0.0214142872 -0.007607369
## 5 2014-10-01 -0.0235734474 -0.0246091992 -0.0243085097 -0.006759468
## 6 2015-01-01 0.0080123694 0.0072758553 0.0074868345 -0.002520833meltedforecasterrors <- melt(forcast_errors, id.vars = 'date')
ggplot(data = meltedforecasterrors, aes(x = date, y = value, color = variable)) + geom_line()
head(outsampleforecast)## date seriesforcast holdoutsample gbm_mle gbm_mm gbm_gmm
## 1 2013-10-01 0.014950458 0.020873837 0.02138044 0.02121181 0.02130124
## 2 2014-01-01 0.041556938 0.046189182 0.02159643 0.02260100 0.02149926
## 3 2014-04-01 0.021486587 0.012073571 0.02182493 0.02187976 0.02192224
## 4 2014-07-01 0.008619872 0.001012503 0.02096821 0.02242679 0.02206379
## 5 2014-10-01 0.004449968 -0.002309500 0.02126395 0.02199901 0.02229970
## 6 2015-01-01 0.032262685 0.029741852 0.02172948 0.02225502 0.02246600# Estimate Root Mean Squared Error (RMSE)
tab_rmse <- data.frame(rbind(GBM.MLE = rmse(outsampleforecast$holdoutsample,outsampleforecast$gbm_mle),
GBM.GMM = rmse(outsampleforecast$holdoutsample,outsampleforecast$gbm_gmm),
GBM.MM = rmse(outsampleforecast$holdoutsample,outsampleforecast$gbm_mm),
ARMA_EGARCH = rmse(outsampleforecast$holdoutsample,outsampleforecast$seriesforcast)))
colnames(tab_rmse) <- c('RMSE')
tab_rmse## RMSE
## GBM.MLE 0.018025225
## GBM.GMM 0.018505961
## GBM.MM 0.019647263
## ARMA_EGARCH 0.006320476# Estimate Mean Absolute Error (MAE)
tab_mae <- data.frame(rbind(GBM.MLE = mae(outsampleforecast$holdoutsample,outsampleforecast$gbm_mle),
GBM.GMM = mae(outsampleforecast$holdoutsample,outsampleforecast$gbm_gmm),
GBM.MM = mae(outsampleforecast$holdoutsample,outsampleforecast$gbm_mm),
ARMA_EGARCH = mae(outsampleforecast$holdoutsample,outsampleforecast$seriesforcast)))
colnames(tab_mae) <- c('MAE')
tab_mae## MAE
## GBM.MLE 0.01624132
## GBM.GMM 0.01667028
## GBM.MM 0.01764616
## ARMA_EGARCH 0.00545459rmse_mae <- cbind(tab_rmse, tab_mae)
rmse_mae## RMSE MAE
## GBM.MLE 0.018025225 0.01624132
## GBM.GMM 0.018505961 0.01667028
## GBM.MM 0.019647263 0.01764616
## ARMA_EGARCH 0.006320476 0.00545459# Diebold & Mariano Test (dmt)
# estimate corresponding p_value
dm_pvalue <- data.frame(rbind(MLE_GMM.pvalue = dm.test(forcast_errors$MLE,
forcast_errors$GMM, alternative = "two.sided", h = 1, power = 2)$p.value,
MLE_MM.pvalue = dm.test(forcast_errors$MLE,forcast_errors$MM, alternative = "two.sided", h = 1, power = 2)$p.value,
MLE_ARMAEGARCH.pvalue = dm.test(forcast_errors$MLE,forcast_errors$ARMA_EGARCH, alternative = "two.sided", h = 1, power = 2)$p.value,
GMM_MM.pvalue = dm.test(forcast_errors$GMM,forcast_errors$MM, alternative = "two.sided", h = 1, power = 2)$p.value,
GMM_ARMAEGARCH.pvalue = dm.test(forcast_errors$GMM,forcast_errors$ARMA_EGARCH, alternative = "two.sided", h = 1, power = 2)$p.value,
MM_ARMAEGARCH.pvalue = dm.test(forcast_errors$MM,forcast_errors$ARMA_EGARCH, alternative = "two.sided", h = 1, power = 2)$p.value))
colnames(dm_pvalue) <- c('P-Value')
results <- list(ArmaEgarch.parameters = model.garch.fit@fit$robust.matcoef, Gbm.parameters = gbm.est, MAE.RMSE = rmse_mae, Diebold.Mariano = dm_pvalue)
# output of analysis
results## $ArmaEgarch.parameters
## Estimate Std. Error t value Pr(>|t|)
## mu -0.002373322 1.529654e-06 -1551.5413278 0.000000e+00
## ar1 -0.091661654 1.822190e-03 -50.3030259 0.000000e+00
## ar2 -0.119877071 8.506420e-04 -140.9254029 0.000000e+00
## ar3 -0.108157906 2.352414e-03 -45.9774163 0.000000e+00
## ar4 0.891410859 7.415993e-03 120.2011514 0.000000e+00
## ma1 0.990996126 1.961266e-04 5052.8388534 0.000000e+00
## ma2 0.952121683 4.419577e-04 2154.3276762 0.000000e+00
## ma3 1.040280467 3.612554e-04 2879.6260846 0.000000e+00
## omega -3.583975692 1.138476e+00 -3.1480476 1.643649e-03
## alpha1 0.076016456 7.675645e-02 0.9903592 3.219986e-01
## beta1 0.574464515 1.382349e-01 4.1557118 3.242764e-05
## gamma1 0.321089259 1.480378e-01 2.1689681 3.008511e-02
##
## $Gbm.parameters
## mu sigma
## MLE 0.07352602 0.05006677
## GMM 0.03724100 0.01500946
## MM 0.07285854 0.04887302
##
## $MAE.RMSE
## RMSE MAE
## GBM.MLE 0.018025225 0.01624132
## GBM.GMM 0.018505961 0.01667028
## GBM.MM 0.019647263 0.01764616
## ARMA_EGARCH 0.006320476 0.00545459
##
## $Diebold.Mariano
## P-Value
## MLE_GMM.pvalue 0.0301716365
## MLE_MM.pvalue 0.0008147499
## MLE_ARMAEGARCH.pvalue 0.0003059662
## GMM_MM.pvalue 0.0093414735
## GMM_ARMAEGARCH.pvalue 0.0001741098
## MM_ARMAEGARCH.pvalue 0.0001825320