Act3&4

a. Background

Briefly describe the selected dependent variable’s background (i.e.,

What is the consumer price index?, What is the consumer sentiment). ## How has been the selected variable’s performance over the period 2005 - 2021?

library(ggplot2)
library(forecast)
library(tseries)
library(zoo)
library(lubridate)
library(TSA)
library(xts)
library(dplyr)
library(scales)
library(vars)
library(readxl)
library(car)
library(glmnet)
library(gridExtra)

ts_data <- read.csv("/Users/pablosancho/Desktop/femsa_ts_data.csv")

str(ts_data)

## 'data.frame':    66 obs. of  10 variables:
##  $ date               : chr  "2005/01" "2005/02" "2005/03" "2005/04" ...
##  $ mexico_cgdp_13     : num  13354788 14104834 13782144 14306524 14107960 ...
##  $ unemployment_rate  : num  0.036 0.035 0.036 0.028 0.033 0.033 0.04 0.033 0.037 0.033 ...
##  $ exchange_rate      : num  11.1 10.8 10.8 10.6 10.7 ...
##  $ mexicp_ipc         : num  12677 13486 16120 17803 19273 ...
##  $ inflation_rate     : num  0.0079 0.008 0.0172 0.0333 0.0087 0.0065 0.0247 0.0405 0.0102 0.0058 ...
##  $ consumer_sentimenta: num  42.2 40.9 41.6 43.6 44.9 ...
##  $ consumer_sentimentb: num  41.7 41.1 41.9 43.6 44.5 ...
##  $ femsa_stock        : num  14.6 16.7 19.5 20 25.9 ...
##  $ igae               : num  83.6 86.4 84.6 89.2 89.5 ...

Change format of date so that it shows quarterly data

ts_data <- ts_data %>%
  mutate(date = sub("/01$", "-01-01", 
               sub("/02$", "-04-01", 
               sub("/03$", "-07-01", 
               sub("/04$", "-10-01", date)))))

ts_data$date <- as.Date(ts_data$date, format = "%Y-%m-%d")

str(ts_data)

## 'data.frame':    66 obs. of  10 variables:
##  $ date               : Date, format: "2005-01-01" "2005-04-01" ...
##  $ mexico_cgdp_13     : num  13354788 14104834 13782144 14306524 14107960 ...
##  $ unemployment_rate  : num  0.036 0.035 0.036 0.028 0.033 0.033 0.04 0.033 0.037 0.033 ...
##  $ exchange_rate      : Time-Series  from 1 to 17.2: 11.1 10.8 10.8 10.6 10.7 ...
##  $ mexicp_ipc         : num  12677 13486 16120 17803 19273 ...
##  $ inflation_rate     : num  0.0079 0.008 0.0172 0.0333 0.0087 0.0065 0.0247 0.0405 0.0102 0.0058 ...
##  $ consumer_sentimenta: num  42.2 40.9 41.6 43.6 44.9 ...
##  $ consumer_sentimentb: num  41.7 41.1 41.9 43.6 44.5 ...
##  $ femsa_stock        : num  14.6 16.7 19.5 20 25.9 ...
##  $ igae               : num  83.6 86.4 84.6 89.2 89.5 ...

b. Visualization

Plot the selected dependent variable (exchange_rate) using a time series format.

# Dependent variable
ggplot(ts_data, aes(x=date, y=exchange_rate)) + 
  geom_line() + 
  labs(title='Exchange Rate Over Time', x='Time', y='Exchange Rate') + 
  theme(axis.text.x = element_text(angle=45, hjust=1))

## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.

Plot the selected explanatory variable (mexicp_ipc) using a time series format

# Exploratory Variable
ggplot(ts_data, aes(x=date, y=mexicp_ipc)) + 
  geom_line() + 
  labs(title='Mexican Stock Exchange Over Time', x='Time', y='IPC') + 
  theme(axis.text.x = element_text(angle=45, hjust=1))

Decompose the selected time series data (both DV and EV) in observed, trend, seasonality, and random.

# Dependent variable exchange_rate
ts_data_ts <- ts(ts_data$exchange_rate, frequency=4, start=c(2005, 1))
decomposed <- decompose(ts_data_ts, type="additive")

plot(decomposed)

### Do the selected time series data show a trend?

• Exchange_rate shows a positive and growing trend.

Do the selected time series data show seasonality?

• The variable shows seasonality.

How is the change of the seasonal component over time?

• The seasonal component has remained stable over time.

# Exploratory variable exchange_rate

ts_data_ts <- ts(ts_data$mexicp_ipc, frequency=4, start=c(2005, 1))
decomposed <- decompose(ts_data_ts, type="additive")

plot(decomposed)

### Do the selected time series data show a trend?

• mexicp_ipc shows a positive and growing trend.

Do the selected time series data show seasonality?

• The variable shows seasonality.

How is the change of the seasonal component over time?

• The seasonal component has remained stable over time.

c. Estimation

Detect if the selected time series data, both DV and EV, are stationary.

Exchange_rate adf test

adf_test <- adf.test(ts_data$exchange_rate, alternative="stationary")

print(paste("ADF Statistic: ", adf_test$statistic))

## [1] "ADF Statistic:  -2.09559986004109"

print(paste("p-value: ", adf_test$p.value))

## [1] "p-value:  0.536267805257164"

The adf test shows a p-value over 5%, meaning the variable is not stationary

exchange_rate_diff1 <- diff(ts_data$exchange_rate)

ts_data$exchange_rate_diff1 <- c(NA, exchange_rate_diff1)

exchange_rate_diff1_clean <- na.omit(ts_data$exchange_rate_diff1)

adf_test_exchange_rate_diff1 <- adf.test(exchange_rate_diff1_clean, alternative = "stationary")

print(paste("ADF Statistic:", adf_test_exchange_rate_diff1$statistic))

## [1] "ADF Statistic: -4.81015236805787"

print(paste("p-value:", adf_test_exchange_rate_diff1$p.value))

## [1] "p-value: 0.01"

print(adf_test_exchange_rate_diff1)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  exchange_rate_diff1_clean
## Dickey-Fuller = -4.8102, Lag order = 3, p-value = 0.01
## alternative hypothesis: stationary

Changed the variable exchange_rate to differenced values to make it stationary.

mexicp_ipc

adf_test_ipc <- adf.test(ts_data$mexicp_ipc, alternative="stationary")

print(paste("ADF Statistic: ", adf_test_ipc$statistic))

## [1] "ADF Statistic:  -2.32023879948318"

print(paste("p-value: ", adf_test_ipc$p.value))

## [1] "p-value:  0.445136389661997"

The adf test shows a p-value over 5%, meaning the variable is not stationary

mexicp_ipc_diff <- diff(ts_data$mexicp_ipc, differences = 1)

ts_data <- ts_data[-1, ] 
ts_data$mexicp_ipc_diff <- mexicp_ipc_diff

# Perform the Augmented Dickey-Fuller test
adf_result <- adf.test(mexicp_ipc_diff, alternative = "stationary")

# Print the result
print(adf_result)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  mexicp_ipc_diff
## Dickey-Fuller = -4.2991, Lag order = 3, p-value = 0.01
## alternative hypothesis: stationary

Changed the variable exchange_rate to differenced values to make it stationary.

Detect if the selected time series data, both DV and EV, show serial autocorrelation.

acf(ts_data$exchange_rate, main='Autocorrelation Plot of Exchange Rate', lag.max=50, ci.col="blue", ci.type="ma")

The variable exchange_rate shows autocorrelation in the first 5 lags. As the lags increase, autocorrelation decreases into the convidence intervals, meaning no autocorrelation.

acf(ts_data$mexicp_ipc, main='Autocorrelation Plot of Mexican IPC', lag.max=50, ci.col="blue", ci.type="ma")

The variable exchange_rate shows autocorrelation in the first 4 lags. As the lags increase, autocorrelation decreases into the convidence intervals, meaning no autocorrelation.

Estimate 3 different time series regression models. You might want to consider ARMA (p,q) and / or ARIMA (p,d,q).

arma_model <- arma(ts_data$exchange_rate, order = c(1, 1))
summary(arma_model)

## 
## Call:
## arma(x = ts_data$exchange_rate, order = c(1, 1))
## 
## Model:
## ARMA(1,1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8096 -0.5677 -0.1481  0.5092  3.1398 
## 
## Coefficient(s):
##            Estimate  Std. Error  t value Pr(>|t|)    
## ar1          0.9815      0.0279   35.181   <2e-16 ***
## ma1         -0.1233      0.1344   -0.918    0.359    
## intercept    0.4221      0.4272    0.988    0.323    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Fit:
## sigma^2 estimated as 0.8054,  Conditional Sum-of-Squares = 50.74,  AIC = 176.39

arima_model <- arima(ts_data$exchange_rate, order = c(1, 1, 1))
summary(arima_model)

## 
## Call:
## arima(x = ts_data$exchange_rate, order = c(1, 1, 1))
## 
## Coefficients:
##           ar1     ma1
##       -0.6456  0.5432
## s.e.   0.3581  0.3804
## 
## sigma^2 estimated as 0.8186:  log likelihood = -84.42,  aic = 172.84
## 
## Training set error measures:
##               ME RMSE MAE MPE MAPE
## Training set NaN  NaN NaN NaN  NaN

arima_model2 <- arima(ts_data$exchange_rate, order = c(2, 1, 2))
summary(arima_model2)

## 
## Call:
## arima(x = ts_data$exchange_rate, order = c(2, 1, 2))
## 
## Coefficients:
##           ar1     ar2      ma1      ma2
##       -0.0243  0.4196  -0.1006  -0.3904
## s.e.   0.6867  0.4326   0.6753   0.3855
## 
## sigma^2 estimated as 0.8167:  log likelihood = -84.34,  aic = 176.69
## 
## Training set error measures:
##               ME RMSE MAE MPE MAPE
## Training set NaN  NaN NaN NaN  NaN

By specifying both DV and EV, estimate a Vector Autoregressive Model (VAR).

var_data <- ts_data[, c("exchange_rate_diff1", "mexicp_ipc_diff")]

# Fitting the VAR model with optimal lag length determined by AIC
var_model <- VAR(var_data, p = VARselect(var_data, lag.max = 10, type = "const")$selection[1], type = "const")

summary(var_model)

## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: exchange_rate_diff1, mexicp_ipc_diff 
## Deterministic variables: const 
## Sample size: 64 
## Log Likelihood: -668.238 
## Roots of the characteristic polynomial:
## 0.06454 0.06454
## Call:
## VAR(y = var_data, p = VARselect(var_data, lag.max = 10, type = "const")$selection[1], 
##     type = "const")
## 
## 
## Estimation results for equation exchange_rate_diff1: 
## ==================================================== 
## exchange_rate_diff1 = exchange_rate_diff1.l1 + mexicp_ipc_diff.l1 + const 
## 
##                          Estimate Std. Error t value Pr(>|t|)  
## exchange_rate_diff1.l1 -2.679e-01  1.528e-01  -1.753   0.0846 .
## mexicp_ipc_diff.l1     -7.766e-05  4.904e-05  -1.584   0.1185  
## const                   2.261e-01  1.205e-01   1.876   0.0654 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.8975 on 61 degrees of freedom
## Multiple R-Squared: 0.05529, Adjusted R-squared: 0.02431 
## F-statistic: 1.785 on 2 and 61 DF,  p-value: 0.1765 
## 
## 
## Estimation results for equation mexicp_ipc_diff: 
## ================================================ 
## mexicp_ipc_diff = exchange_rate_diff1.l1 + mexicp_ipc_diff.l1 + const 
## 
##                        Estimate Std. Error t value Pr(>|t|)
## exchange_rate_diff1.l1 560.4056   484.4991   1.157    0.252
## mexicp_ipc_diff.l1       0.1469     0.1555   0.945    0.348
## const                  411.5259   381.9499   1.077    0.286
## 
## 
## Residual standard error: 2846 on 61 degrees of freedom
## Multiple R-Squared: 0.0233,  Adjusted R-squared: -0.008721 
## F-statistic: 0.7277 on 2 and 61 DF,  p-value: 0.4872 
## 
## 
## 
## Covariance matrix of residuals:
##                     exchange_rate_diff1 mexicp_ipc_diff
## exchange_rate_diff1              0.8054           -1448
## mexicp_ipc_diff              -1448.0533         8096965
## 
## Correlation matrix of residuals:
##                     exchange_rate_diff1 mexicp_ipc_diff
## exchange_rate_diff1               1.000          -0.567
## mexicp_ipc_diff                  -0.567           1.000

# Calculate the AIC of the fitted VAR model
aic_value <- AIC(var_model)

# Print the AIC value
print(aic_value)

## [1] 1348.476

Estimate a Granger Causality Test. Briefly, interpret the estimated results.

# Granger causality test for "exchange_rate" causing "mexicp_ipc"
granger_test1 <- causality(var_model, cause = "exchange_rate_diff1")
print(granger_test1)

## $Granger
## 
##  Granger causality H0: exchange_rate_diff1 do not Granger-cause
##  mexicp_ipc_diff
## 
## data:  VAR object var_model
## F-Test = 1.3379, df1 = 1, df2 = 122, p-value = 0.2497
## 
## 
## $Instant
## 
##  H0: No instantaneous causality between: exchange_rate_diff1 and
##  mexicp_ipc_diff
## 
## data:  VAR object var_model
## Chi-squared = 15.571, df = 1, p-value = 7.945e-05

granger_test2 <- causality(var_model, cause = "mexicp_ipc_diff")
print(granger_test2)

## $Granger
## 
##  Granger causality H0: mexicp_ipc_diff do not Granger-cause
##  exchange_rate_diff1
## 
## data:  VAR object var_model
## F-Test = 2.5076, df1 = 1, df2 = 122, p-value = 0.1159
## 
## 
## $Instant
## 
##  H0: No instantaneous causality between: mexicp_ipc_diff and
##  exchange_rate_diff1
## 
## data:  VAR object var_model
## Chi-squared = 15.571, df = 1, p-value = 7.945e-05

Exchange rate does not granger-cause mexicp_ipc and viceversa. There is no existing causality meaning we fail to reject the null hypothesis.

d. Evaluation

Based on diagnostic tests and RMSE and / or AIC, compare the 3 estimated time series regression models from Assignment 3 and the estimated VAR model results, and select the time series model that better fits the data to estimate the forecast.

# Creating a data frame with the model names and their corresponding AIC values
aic_comparison <- data.frame(
  Model = c("ARMA(1,1)", "ARIMA(1,1,1)", "ARIMA2(2,1,2)", "VAR"),
  AIC = c(176.39, 172.84, 176.69, 1348.476)
)

# Displaying the table
print(aic_comparison)

##           Model      AIC
## 1     ARMA(1,1)  176.390
## 2  ARIMA(1,1,1)  172.840
## 3 ARIMA2(2,1,2)  176.690
## 4           VAR 1348.476

The selected model is Arima (1, 1, 1) because it has the lowest AIC, with 172.84.

e. Forecast

By using the selected model, make a forecast for the next 6 periods. In doing so, include a time series plot showing your forecast.

Arima_forecast <-forecast(ts_data$exchange_rate, h=6)
plot(Arima_forecast)

autoplot(Arima_forecast)