<jupyter_start><jupyter_text>Exercise - Dynamic models and model selection**(1)** Fit the following ARDL models:> **(a)** $y_t = \alpha + \beta_1 y_{t-1} + \beta_2 ipr_t + \beta_2 ipr_{t-1} + \beta_3 su_t + \beta_4 su_{t-1} + \beta_5 sr_t + \beta_6 sr_{t-1} + \epsilon_T$> **(b)** $y_t = \alpha + \beta_1 y_{t-1} + \epsilon_t$> **(c)** $y_t = \alpha + \beta_1 ydiff_{t-1} + \beta_2 iprdiff_t + \beta_3 sudiff_t + \beta_4 srdiff_t + \epsilon_t$> **(d)** $y_t = \alpha + \beta_1 ipr_{t-1} + \beta_2 su_{t-1} + \beta_3 sr_{t-1} + \epsilon_t$> **(e)** $y_t = \alpha + \beta_1 ipr_T + \beta_2 ipr_{t-1} + \beta_3 su_t + \beta_4 su_{t-1} + \beta_5 sr_t + \beta_6 sr_{t-1} + \epsilon_T$> **(f)** $y_t = \alpha + \beta_1 y_{t-1} + \beta_2 ipr_t + \beta_3 su_t + \beta_4 sr_t + \epsilon_t$> **(g)** $y_t = \alpha + \beta_1 ipr_{t-1} + \beta_2 su_{t-1} + \beta_3 sr_{t-1} + \epsilon_t$> **(h)** $y_t = \alpha + \beta_1 ipr_t + \beta_2 su_t + \beta_3 sr_t + \epsilon_t$ onde $\epsilon_t = \rho \epsilon_{t-1} + u_{t-1}$**(2)** Test the optimal ARDL model specification chosen before in terms of:> **(a)** AIC and BIC> **(b)** Serial autocorrelation > **(c)** Homoskedasticity**(3)** Use the optimal model to make:> **(a)** Static forecast > **(b)** Recursive forecast> **(c)** Evaluate the forecasts in terms of RMSE and MAE<jupyter_code>library(urca)
library(readxl)
library(vars)
library(tsDyn)
library(nlme)
library(data.table)
library(ggplot2)
library(lmtest)
library(tseries)
library(xtable)
library(forecast)
library(Metrics)
library(here)
library(reshape2)
library(Metrics)
us_gdp = fread(here('src', 'data', 'ex2_dynamic_gdp_us.csv'))
us_gdp$d1ipr=c(NA,diff(us_gdp$ipr,differences=1)[1:length(us_gdp$ipr)-1])
us_gdp$d1su=c(NA,diff(us_gdp$su,differences=1)[1:length(us_gdp$su)-1])
us_gdp$d1sr=c(NA,diff(us_gdp$sr,differences=1)[1:length(us_gdp$sr)-1])
head(us_gdp)
dim(us_gdp)<jupyter_output><empty_output><jupyter_text>**(1)** Fit the following ARDL models:<jupyter_code>model_list = list()
aic_bic_list = list()
l = list()
model_list$ardl1111 = lm(y ~ y1 + ipr + ipr1 + su + su1 + sr + sr1, us_gdp[2:70] )

l$AIC = AIC(model_list$ardl1111)/model_list$ardl1111$df.residual
l$BIC = BIC(model_list$ardl1111)/model_list$ardl1111$df.residual

aic_bic_list[['arfl1111']] = l
l = list()
model_list$ar1 = lm(y ~ y1 , us_gdp[2:70] )

l$AIC = AIC(model_list$ar1)/model_list$ar1$df.residual
l$BIC = BIC(model_list$ar1)/model_list$ar1$df.residual

aic_bic_list[['ar1']] = l
l = list()
model_list$diffmodel = lm(y ~ offset(1*y1)+ d1ipr + d1su + d1sr, data = us_gdp[2:70])

l$AIC = AIC(model_list$diffmodel)/model_list$diffmodel$df.residual
l$BIC = BIC(model_list$diffmodel)/model_list$diffmodel$df.residual

aic_bic_list[['diffmodel']] = l
l = list()
model_list$leadingindcmodels = lm(y ~ ipr1 + su1 + sr1, data = us_gdp[2:70])

l$AIC = AIC(model_list$leadingindcmodels)/model_list$leadingindcmodels$df.residual
l$BIC = BIC(model_list$leadingindcmodels)/model_list$leadingindcmodels$df.residual

aic_bic_list[['leadingindcmodels']] = l
l = list()
model_list$distlagsmodel = lm(y ~ ipr + ipr1 + su + su1 + sr + sr1, data = us_gdp[2:70])

l$AIC = AIC(model_list$distlagsmodel)/model_list$distlagsmodel$df.residual
l$BIC = BIC(model_list$distlagsmodel)/model_list$distlagsmodel$df.residual

aic_bic_list[['distlagsmodel']] = l
l = list()
model_list$partialadjmodel = lm(y ~ y1 + ipr + su + sr, data = us_gdp[2:70])

l$AIC = AIC(model_list$partialadjmodel)/model_list$partialadjmodel$df.residual
l$BIC = BIC(model_list$partialadjmodel)/model_list$partialadjmodel$df.residual

aic_bic_list[['partialadjmodel']] = l
l = list()
model_list$deadstartmodel = lm(y ~ y1 + ipr1+su1+sr1, data = us_gdp[2:70])

l$AIC = AIC(model_list$deadstartmodel)/model_list$deadstartmodel$df.residual
l$BIC = BIC(model_list$deadstartmodel)/model_list$deadstartmodel$df.residual

aic_bic_list[['deadstartmodel']] = l
l = list()
model_list$staticar1errormodel = gls(y ~ ipr+ su + sr, correlation =corARMA(p=1,q=0), data = us_gdp[2:70])


l$AIC = AIC(model_list$staticar1errormodel)/model_list$staticar1errormodel$dims$N
l$BIC = BIC(model_list$staticar1errormodel)/model_list$staticar1errormodel$dims$N

aic_bic_list[['staticar1errormodel']] = l<jupyter_output><empty_output><jupyter_text>**(2)** Test the optimal ARDL model specification chosen before in terms of: **(a)** AIC and BIC<jupyter_code>do.call('cbind', aic_bic_list)<jupyter_output><empty_output><jupyter_text>Its possible to note that both in terms of AIC and BIC, the "partialadjmodel" is the best fit between the available models.<jupyter_code>name = 'ardl1111'
modd = model_list[[name]]
test = bgtest(modd)
test$p.value<jupyter_output><empty_output><jupyter_text>**(b)** Serial autocorrelation  Breusch-Godfrey test<jupyter_code>bg_pval_list = list()
for (name in names(model_list)){
    if (name != 'staticar1errormodel'){
        modd = model_list[[name]]
        test = bgtest(modd)
        bg_pval_list[[name]] = test$p.value
    }
}
do.call('cbind', bg_pval_list)<jupyter_output><empty_output><jupyter_text>**(c)** Homoskedasticity Breusch-Pagan-Godfrey test<jupyter_code>bp_pval_list = list()
for (name in names(model_list)){
    if (name != 'staticar1errormodel'){
        modd = model_list[[name]]
        test = bptest(modd)
        bp_pval_list[[name]] = test$p.value
    }
}
do.call('cbind', bg_pval_list)<jupyter_output><empty_output><jupyter_text>**(3)** Use the optimal model to make: **(a)** Static forecast<jupyter_code>yhat_list = list()
yse_list = list()
ylow_list = list()
yup_list = list()
yhat_static_list = list()
yse_static_list = list()
ylow_static_list = list()
yup_static_list = list()
for (name in names(model_list)){    
    modd = model_list[[name]]
    yhat = predict(modd, newdata = us_gdp[71:114])
    yse = sqrt(sum(modd$residuals^2) / 194)
    yup = yhat + 1.96 * yse
    ylow = yhat - 1.96 * yse
    
    yhat_static_list[[name]] = yhat
    yse_static_list[[name]] = yse
    ylow_static_list[[name]] = ylow
    yup_static_list[[name]] = yup
}
yhat_list$static = yhat_static_list
yse_list$static = yse_static_list
ylow_list$static = ylow_static_list
yup_list$static = yup_static_list<jupyter_output><empty_output><jupyter_text>**(b)** Recursive forecast<jupyter_code>yhat_rec_list = list()
yse_rec_list = list()
ylow_rec_list = list()
yup_rec_list = list()

for (name in names(model_list)){
    yrec = rep(0, 43)
    yrec_se = rep(0, 43)
    yrec_up = rep(0, 43)
    yrec_low = rep(0, 43)
    j=0
    for(i in 70:113) {
        if (name == 'ardl1111'){
            modd = lm(y ~ y1 + ipr + ipr1 + su + su1 + sr + sr1, us_gdp[1:i] )
        }
        else if (name == 'ar1'){
            modd = lm(y ~ y1 , us_gdp[2:70] )
        }
        else if (name == 'diffmodel'){
            modd = lm(y ~ offset(1*y1)+ d1ipr + d1su + d1sr, data = us_gdp[1:i])
        } 
        else if (name == 'leadingindcmodels'){
            modd = lm(y ~ ipr1 + su1 + sr1, data = us_gdp[1:i])
        }
        else if (name == 'distlagsmodel'){
            modd = lm(y ~ ipr + ipr1 + su + su1 + sr + sr1, data = us_gdp[1:i])
        }
        else if (name == 'partialadjmodel'){
            modd = lm(y ~ y1 + ipr + su + sr, data = us_gdp[1:i])
        }
        else if (name == 'deadstartmodel'){
            modd = lm(y ~ y1 + ipr1+su1+sr1, data = us_gdp[1:i])
        }else{
            modd = gls(y ~ ipr+ su + sr, correlation =corARMA(p=1,q=0), data = us_gdp[1:i])
        }
        
        yrec[j] = predict(modd, newdata = us_gdp[(i+1)])
        yrec_se[j] = sqrt(sum(modd$residuals^2 / 194))
        yrec_up[j] = yrec[j] + 1.96 * yrec_se[j]
        yrec_low[j] = yrec[j] - 1.96 * yrec_se[j]
        j = j + 1
    }
    yhat_rec_list[[name]] = yrec
    yse_rec_list[[name]] = yrec_se
    ylow_rec_list[[name]] = yrec_up
    yup_rec_list[[name]] = yrec_low
}

yhat_list$rec = yhat_rec_list
yse_list$rec = yse_rec_list
ylow_list$rec = ylow_rec_list
yup_list$rec = yup_rec_list
df_pred_static = as.data.frame(cbind(x=seq(0,43), y1=as.numeric(us_gdp$y[71:114]), apply(as.data.frame(do.call('cbind', yhat_list$static)), 2, function(x) as.numeric(x))))
head(df_pred_static)
df_pred_static_melt = melt(df_pred_static, id = 'x')
ggplot(data=df_pred_static_melt, aes(x=x, y=value, colour=variable)) + geom_line(size = 0.5)
df_pred_rec = as.data.frame(cbind(x=seq(0,43), y1=us_gdp$y[71:114], do.call('cbind', yhat_list$rec)))
head(df_pred_rec)
df_pred_rec_melt = melt(df_pred_rec, id = "x")
ggplot(data=df_pred_rec_melt, aes(x=x, y=value, colour=variable)) + geom_line(size = 0.5)<jupyter_output><empty_output><jupyter_text>**(c)** Evaluate the forecasts in terms of RMSE and MAE<jupyter_code>static_list = list()

mae_list = list()
rmse_list = list()
for (col in colnames(df_pred_static)){
    if ((col != 'x') & (col != 'y1')){
        mae_list[[col]] = mae(actual = df_pred_static$y1, predicted = df_pred_static[[col]])
        rmse_list[[col]] = rmse(actual = df_pred_static$y1, predicted = df_pred_static[[col]])
        
        static_list$mae = mae_list
        static_list$rmse = rmse_list
    }
}
do.call('cbind', static_list)
rec_list = list()

mae_list = list()
rmse_list = list()
for (col in colnames(df_pred_rec)){
    if ((col != 'x') & (col != 'y1')){
        mae_list[[col]] = mae(actual = df_pred_rec$y1, predicted = df_pred_rec[[col]])
        rmse_list[[col]] = rmse(actual = df_pred_rec$y1, predicted = df_pred_rec[[col]])
        
        rec_list$mae = mae_list
        rec_list$rmse = rmse_list
    }
}
do.call('cbind', rec_list)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1667016806722689, 'fraction_numerical': 0.028466386554621848, 'mean_word_length': 4.620425029515939, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 51, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '528075', 'license': 'MIT', 'n_tokens_mistral': 4243, 'n_tokens_neox': 4072, 'n_words': 710, 'path': 'src/dynamic_model_n_model_selection.ipynb', 'repo_name': 'dcuoliveira/FGV_Forecasting', 'size': '9520'}