<jupyter_start><jupyter_text>Exercise 1**(1)** Forecasting with linear models:> **(a)** Estimate four linear models unsing the OLS estimator> **(b)** Forecast n steps ahead using the estimated models> **(c)** Forecast n steps ahead (recursively) using the estimated models> **(d)** Compute confidence intervals for both **(b)** and **(c)** forecasts<jupyter_code>library(data.table)
library(readr)
library(here)<jupyter_output><empty_output><jupyter_text>(a) Estimate four linear models unsing the OLS estimator The linear regression model (LRM) can be written as follows:$\mathbf{y} = X\mathbf{\beta} + \mathbf{\epsilon}$where $\mathbf{y}$ is a Tx1 vector of dependent variables, $X$ is a Txk matrix of independen variables, $\mathbf{\beta}$ is a kx1 vector of parameters, and $\mathbf{\epsilon}$ is a Tx1 vector of independent error terms.The LRM assumptions are:> A1. $\mathbf{E}[\mathbf{\epsilon}]=0$> A2. $\mathbf{E}[\mathbf{\epsilon}^{'}\mathbf{\epsilon}]=\sigma^2\mathbf{I}$> A3. $X \perp \!\!\! \perp \mathbf{\epsilon}$> A4. $X^{'}X$ is non singular> A5. $X$ is weekly stationaryFurthermore, the parameters $\mathbf{\beta}$ can be estimated using the ordinary least squares (OLS) method. This method seeks to find the vector $\mathbf{\beta}$ that solves the following optimization problem:$argmax_{\mathbf{\beta}} SSE(X, \mathbf{y}; \mathbf{\beta}) \quad SSE = \mathbf{\epsilon}^{'}\mathbf{\epsilon} = (\mathbf{y} - X\mathbf{\beta})^{'}(\mathbf{y} - X\mathbf{\beta})$which has the following solution:$\hat{\mathbf{\beta}} = (X^{'}X)^{-1}X^{'}\mathbf{y}$For the variance, the OLS estimator ir given by:$\hat{\mathbf{\sigma}}^2 = \frac{\hat{\mathbf{\epsilon}}^{'}\hat{\mathbf{\epsilon}}}{T-k} \quad where \quad \hat{\mathbf{\epsilon}} = \mathbf{y} - X\hat{\mathbf{\beta}}$It is possible to show that, under the assumtions A1-A5, the OLS estimator is the best linear unbiased estimator and has the lowest variance in this class.<jupyter_code>df = readr::read_csv(here("src", "data", "ex2_regress_gdp.csv"))
head(df)
par(mfrow=c(3, 2))
plot(df$y, type = 'line', main = 'GDP')
plot(df$pr, type = 'line', main = 'Industrial production')
plot(df$sr, type = 'line', main = 'Consumer sentiment index')
plot(df$su, type = 'line', main = 'Eurostoxx')
plot(df$ipr, type = 'line', main = 'inflation rate')<jupyter_output>Warning message in plot.xy(xy, type, ...):
"plot type 'line' will be truncated to first character"Warning message in plot.xy(xy, type, ...):
"plot type 'line' will be truncated to first character"Warning message in plot.xy(xy, type, ...):
"plot type 'line' will be truncated to first character"Warning message in plot.xy(xy, type, ...):
"plot type 'line' will be truncated to first character"Warning message in plot.xy(xy, type, ...):
"plot type 'line' will be truncated to first character"<jupyter_text>By doing a quick check on the time series under analysis it seems that some of the LRM hypothesis may be violated. For instance, all the series seems to have a variance that is not constant over time, violating the A2 assumption. Furthermore, the GDP series, in particular, has an outlier in its history, which may undermine the standard error estimates of the OLS.Despite these facts, we estimate four specifications of the LRM as follows:$gdp_t = \beta_1 ipr_t + \beta_2 su_t + \beta_3 pr_t + \beta_4 sr_t + \epsilon_t$$gdp_t = \beta_1 ipr_t + \beta_2 su_t + \beta_3 sr_t + \epsilon_t$$gdp_t = \beta_1 ipr_t + \beta_2 su_t + \epsilon_t$$gdp_t = \beta_1 ipr_t +  \beta_2 pr_t + \beta_3 sr_t + \epsilon_t$where:> $ipr_t$ is the industrial production for the euro area at time $t$> $su_t$ is the Eurostoxx returns at $t$> $pr_t$ is the euro area inflation rate at $t$> $sr_t$ is the euro area consumer sentiment index at $t$All series are in monthly frequency.<jupyter_code>gdp.fit <- list()
gdp.formula <- c('y ~ ipr + su + pr + sr', 'y ~ ipr + su + sr',
                 'y ~ ipr + su', 'y ~ ipr + pr + sr')

for (model in 1:4) {
  gdp.fit[[model]] <- lm(gdp.formula[model], data = df)
}
est_df = df[1:44,] # estimation window
fore_df = df[45:70,] # forecast window

gdp.est <- list()
for (model in 1:4) {
  gdp.est[[model]] <- lm(gdp.formula[model], data = est_df)
  summary(gdp.est[[model]])
}<jupyter_output><empty_output><jupyter_text>The model estimation summaries are:<jupyter_code>summary(gdp.fit[[1]])
summary(gdp.fit[[2]])
summary(gdp.fit[[3]])
summary(gdp.fit[[4]])<jupyter_output><empty_output><jupyter_text>One common statistics for model evaluation is the coefficient of determination or $R^2$:$R^2 = 1-\frac{\mathbf{\epsilon}^{'}\mathbf{\epsilon}}{\mathbf{y}^{'}\mathbf{y}}$The $R^2$ compares the variability in the dependent variable $\mathbf{y}$ explained by the model with the variability of $\mathbf{\epsilon}$, which is the unexplanied amount of variability.Despite been naive to analyse the output of a regression in terms of the $R^2$ only, we see that the first model has the highest between the tested models ($R^2$=0.8008).One clear undesirable property of the $R^2$ is that its an increasing function of the number of parameters, therefore the adjusted $R^2$ is typically a better measure.$\bar{R}^2 = 1-\frac{\mathbf{\epsilon}^{'}\mathbf{\epsilon}}{T-k}\frac{T-1}{\mathbf{y}^{'}\mathbf{y}}$We can note that despite the first model having the highest $R^2$=0.8008, its adjusted version has $\bar{R}^2=0.7886$ which is pretty close to the second model $\bar{R}^2=0.7866$. Therefore, it seems that the unadjusted r-quared of the first model has been inflated by the number of variables and we should the second model instead.  (b) Forecast n steps ahead using the estimated models + (c) Forecast n steps ahead (recursively) using the estimated models The OLS regression forecast is given by:$\hat{\mathbf{y}}_{T+h} = X_{T+h}\hat{\mathbf{\beta}}$and the forecast error is simply:$e_{t+h} = \mathbf{y}_{T+h} - \hat{\mathbf{y}}_{T+h}$<jupyter_code>gdp.fore = list()
gdp.rec = list()

for (model in 1:4) {
  gdp.fore[[model]] = predict(gdp.est[[model]], newdata = fore_df)
  
  gdp.rec[[model]] = rep(0, 26)
  for (i in 1:26) {
    ols.rec = lm(gdp.formula[model], data = df[1:(43 + i),])
    gdp.rec[[model]][i] = predict(ols.rec, newdata = df[44 + i,])
  }
}
plot.label <- 2007:2013
par(mfrow=c(2, 2))
for (model in 1:4) {
  gdp.plot = cbind(data.table(fore_df$y),
                    data.table(gdp.rec[[model]]),
                    data.table(gdp.fore[[model]]))
  setnames(gdp.plot, c('Y', paste0('YRFOREG', model),
                       paste0('YFOREG', model)))
  
  plot(gdp.plot[[1]], type = 'n', xlab = '', ylab = '',
       xaxt = 'n', ylim = c(-3, 2))
  axis(1, at = c(0:6) * 4 + 2, labels = plot.label)
  for (l in 1:3) {
    lines(gdp.plot[[l]], lty = l)
  }
  legend('bottomright', legend = colnames(gdp.plot), inset = 0.05, lty = 1:3)
}
error_fore = list()
error_rec = list()
for (model in 1:4) {
    gdp.plot = cbind(data.table(fore_df$y),
                    data.table(gdp.rec[[model]]),
                    data.table(gdp.fore[[model]]))
    error_fore[model] = sum((gdp.plot[,1]-gdp.plot[,3])^2 / dim(gdp.plot[,3])[1]
    error_rec[model] = sum((gdp.plot[,1]-gdp.plot[,2])^2) / dim(gdp.plot[,2])[1]
}
error_fore
sum(unlist(error_fore))
error_rec
sum(unlist(error_rec))<jupyter_output><empty_output><jupyter_text>It can be seem from the above that model two has the lowest forecast error (0.14307) between the batch forecast method, and the lowest error (0.12327) between the recursive forecast method.Furthermore, it can be seem that, as expected, the recursive method has consistently beaten the batch forecast method in terms of mean squared error across all models. **(d) Compute confidence intervals for both (b) and (c) forecasts** A typical parametric assumption of the LRM is:$\epsilon \sim N(0, \sigma^2 \mathbf{I})$this implies that the distribution of the OLS estimator can be written as:$\sqrt{T}(\hat{\beta}-\beta) \sim N(0, \sigma^2(\frac{X^{'}X}{T})^{-1})$and the forecast errors have the following distribution:$\frac{\mathbf{y}_{T+h} - \hat{\mathbf{y}}_{T+h}}{\sqrt{Var(e_{T+h})}} \sim N(0,1) \implies \mathbf{y}_{T+h} \sim N(\hat{\mathbf{y}}_{T+h}, Var(e_{T+h}))$This means that we can use the above density to construct confidence intervals for our forecasts. In particular, a $[1-\alpha]$% forecast interval is represented by:$[\hat{\mathbf{y}}_{T+h}-c_{\alpha/2}\sqrt{Var(e_{T+h})} ; \hat{\mathbf{y}}_{T+h}+c_{\alpha/2}\sqrt{Var(e_{T+h})}]$ Therefore, for a 95% area interval we have that:<jupyter_code>gdp.fore.ic.se = list()
gdp.fore.ic.up = list()
gdp.fore.ic.low = list()
for (i in 1:length(gdp.est)){
    gdp.fore.ic.se[[i]] = sqrt(sum(gdp.est[[1]]$residuals^2) / dim(fore_df)[1])
    gdp.fore.ic.up[[i]] = gdp.fore[[i]] + 1.96 * unlist(gdp.fore.ic.se)
    gdp.fore.ic.low[[i]] = gdp.fore[[i]] - 1.96 * unlist(gdp.fore.ic.se)
}
par(mfrow=c(2,2))

for (i in 1:length(gdp.est)){
    dat = cbind(gdp.fore.ic.up[[i]], gdp.fore.ic.low[[i]], gdp.fore[[i]], fore_df$y)
    colnames(dat) = names = c('up', 'lower', 'predict', 'real')    
    
    matplot(dat, type = c("b"),pch=1,col = 1:4)
    legend("topleft", legend=names, col=1:4, pch=1)
}<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1730391085864103, 'fraction_numerical': 0.019445051343674895, 'mean_word_length': 4.572124163116251, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 37, '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': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '17729', 'license': 'MIT', 'n_tokens_mistral': 3304, 'n_tokens_neox': 3132, 'n_words': 966, 'path': 'src/Linear_model_confidence_interval.ipynb', 'repo_name': 'dcuoliveira/FGV_Forecasting', 'size': '9154'}