Evidence 2

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)

I. Introduction

Time series analysis is a statistical technique that involves analyzing data points collected or recorded at specific time intervals. This type of analysis is used to identify patterns, trends, cycles, and seasonal variations in the data over time. The primary goal of time series analysis is to forecast future values based on historical data, making it particularly useful in fields such as finance, economics, weather forecasting, and many other areas where data evolves over time (Chatfield, 2003).

Stationarity is a fundamental concept in time series analysis, which implies that the statistical properties of the series, such as mean and variance, remain constant over time. Stationary time series are easier to model and predict because their behavior is consistent over time. Non-stationary data, on the other hand, may exhibit trends, seasonality, or varying variance, complicating the analysis and potentially leading to misleading results. Ensuring stationarity often involves transforming or differencing the data to remove trends and seasonality, allowing for more accurate modeling and forecasting (Box et al., 2015).

Time series analysis is crucial for business intelligence because it enables organizations to make data-driven decisions by forecasting future trends, demand, or performance based on historical data. It can be applied to various business problems, such as predicting sales, stock prices, or customer behavior. This analysis helps businesses anticipate changes, manage inventory, optimize marketing strategies, and improve operational efficiency by understanding patterns over time and responding proactively.

Sources Chatfield, C. (2003). The Analysis of Time Series: An Introduction. Chapman and Hall/CRC. Box, G.E.P., Jenkins, G.M., Reinsel, G.C., & Ljung, G.M. (2015). Time Series Analysis: Forecasting and Control. Wiley.

II. Description of the Problem Situation

For this problem-situation, we are looking to develop a series of models that are effectively able of predicting Coca-Cola’s sales through the Metropolitan Area of Guadalajara, Jalisco. To do this, we are identifying variables that have strong correlation with our dependent variables (sales), in hopes of detecting patterns that will help Coca-Cola’s sales in the future.

III. Data and Methodology

Format transformation

str(cocacola)

## tibble [48 × 15] (S3: tbl_df/tbl/data.frame)
##  $ tperiod           : POSIXct[1:48], format: "2021-01-15" "2021-02-15" ...
##  $ sales_unitboxes   : num [1:48] 5516689 5387496 5886747 6389182 6448275 ...
##  $ consumer_sentiment: num [1:48] 38.1 37.5 38.5 37.8 38 ...
##  $ CPI               : num [1:48] 87.1 87.3 87.6 87.4 87 ...
##  $ inflation_rate    : num [1:48] -0.09 0.19 0.41 -0.26 -0.5 0.17 0.15 0.21 0.37 0.51 ...
##  $ unemp_rate        : num [1:48] 0.0523 0.0531 0.0461 0.051 0.0552 ...
##  $ gdp_percapita     : num [1:48] 11660 11660 11660 11626 11626 ...
##  $ itaee             : num [1:48] 104 104 104 108 108 ...
##  $ itaee_growth      : num [1:48] 0.0497 0.0497 0.0497 0.0318 0.0318 ...
##  $ pop_density       : num [1:48] 98.5 98.5 98.5 98.8 98.8 ...
##  $ job_density       : num [1:48] 18.3 18.5 18.6 18.7 18.7 ...
##  $ pop_minwage       : num [1:48] 9.66 9.66 9.66 9.59 9.59 ...
##  $ exchange_rate     : num [1:48] 14.7 14.9 15.2 15.2 15.3 ...
##  $ max_temperature   : num [1:48] 28 31 29 32 34 32 29 29 29 29 ...
##  $ holiday_month     : num [1:48] 0 0 0 1 0 0 0 0 1 0 ...

cocacola$holiday_month <- as.factor(cocacola$holiday_month)

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

day_month <- format(cocacola$tperiod, "%d-%m")

corrected_year <- as.numeric(substr(day_month, 1, 2)) + 2000
corrected_month_day <- paste0(corrected_year, "-", substr(day_month, 4, 5))

cocacola$tperiod <- as.Date(paste0(corrected_month_day, "-01"), format = "%Y-%m-%d")

head(cocacola$tperiod)

## [1] "2015-01-01" "2015-02-01" "2015-03-01" "2015-04-01" "2015-05-01"
## [6] "2015-06-01"

str(cocacola)

## tibble [48 × 15] (S3: tbl_df/tbl/data.frame)
##  $ tperiod           : Date[1:48], format: "2015-01-01" "2015-02-01" ...
##  $ sales_unitboxes   : num [1:48] 5516689 5387496 5886747 6389182 6448275 ...
##  $ consumer_sentiment: num [1:48] 38.1 37.5 38.5 37.8 38 ...
##  $ CPI               : num [1:48] 87.1 87.3 87.6 87.4 87 ...
##  $ inflation_rate    : num [1:48] -0.09 0.19 0.41 -0.26 -0.5 0.17 0.15 0.21 0.37 0.51 ...
##  $ unemp_rate        : num [1:48] 0.0523 0.0531 0.0461 0.051 0.0552 ...
##  $ gdp_percapita     : num [1:48] 11660 11660 11660 11626 11626 ...
##  $ itaee             : num [1:48] 104 104 104 108 108 ...
##  $ itaee_growth      : num [1:48] 0.0497 0.0497 0.0497 0.0318 0.0318 ...
##  $ pop_density       : num [1:48] 98.5 98.5 98.5 98.8 98.8 ...
##  $ job_density       : num [1:48] 18.3 18.5 18.6 18.7 18.7 ...
##  $ pop_minwage       : num [1:48] 9.66 9.66 9.66 9.59 9.59 ...
##  $ exchange_rate     : num [1:48] 14.7 14.9 15.2 15.2 15.3 ...
##  $ max_temperature   : num [1:48] 28 31 29 32 34 32 29 29 29 29 ...
##  $ holiday_month     : Factor w/ 2 levels "0","1": 1 1 1 2 1 1 1 1 2 1 ...

Dependent variable analysis

ggplot(cocacola, aes(x=tperiod, y=sales_unitboxes)) + 
  geom_line() + 
  labs(title='Box Sales Over Time', x='Time', y='sales_unitboxes') + 
  theme(axis.text.x = element_text(angle=45, hjust=1))

cocacola$Moving_Avg <- rollmean(cocacola$sales_unitboxes, 12, fill = NA)

ggplot(cocacola, aes(x=tperiod)) + 
  geom_line(aes(y=sales_unitboxes, colour="Original")) + 
  geom_line(aes(y=Moving_Avg, colour="4-Periods Moving Average")) +
  labs(title='Moving Average Plot', x='Date', y='Sales (Unit Boxes)') +
  scale_colour_manual("", breaks = c("Original", "12-Periods Moving Average"), values = c("blue", "red"))

cocacola_ts <- ts(cocacola$sales_unitboxes, frequency=12, start=c(2015, 1))
decomposed <- decompose(cocacola_ts, type="additive")

plot(decomposed)

adf_test <- adf.test(cocacola$sales_unitboxes, alternative="stationary")

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

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

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

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

sales_ts <- ts(cocacola$sales_unitboxes, frequency=12) 

acf(sales_ts, main="Autocorrelation of Sales (Unit Boxes)")

pacf(sales_ts, main="Partial Autocorrelation of Sales (Unit Boxes)")

The dependent variable “sales_unitboxes” has a trendline that moves up, which indicates increasing over time. It also displays seasonality, which is likely due to temperature. Finally, we can also conclude it’s stationary, due to the p-value (0.01) and the autocorrelation demonstrated in the ACF plot.

IV. Time Series Regression Analysis

Model Comparison

p_values <- 0:5  
q_values <- 0:5 

best_aic_arma <- Inf
best_order_arma <- c(0, 0)
best_model_arma <- NULL

for (p in p_values) {
  for (q in q_values) {
    model <- tryCatch({
      arima(sales_ts, order = c(p, 0, q))  
    }, error = function(e) NULL)  
    
    if (!is.null(model)) {
      current_aic <- AIC(model)  
      
      if (current_aic < best_aic_arma) {
        best_aic_arma <- current_aic
        best_order_arma <- c(p, q)
        best_model_arma <- model
      }
    }
  }
}

p_values <- 0:5  
d_values <- 1:2  
q_values <- 0:5 

best_aic_arima <- Inf
best_order_arima <- c(0, 1, 0)
best_model_arima <- NULL

for (p in p_values) {
  for (d in d_values) {
    for (q in q_values) {
      model <- tryCatch({
        arima(sales_ts, order = c(p, d, q))  
      }, error = function(e) NULL)  
    
      if (!is.null(model)) {
        current_aic <- AIC(model)  
        
        if (current_aic < best_aic_arima) {
          best_aic_arima <- current_aic
          best_order_arima <- c(p, d, q)
          best_model_arima <- model
        }
      }
    }
  }
}

p_values <- 0:5  

best_aic_ar <- Inf
best_order_ar <- c(0, 0)
best_model_ar <- NULL

for (p in p_values) {
  model <- tryCatch({
    arima(sales_ts, order = c(p, 0, 0))  
  }, error = function(e) NULL)  
    
  if (!is.null(model)) {
    current_aic <- AIC(model)  
    
    if (current_aic < best_aic_ar) {
      best_aic_ar <- current_aic
      best_order_ar <- c(p, 0)
      best_model_ar <- model
    }
  }
}

q_values <- 0:5  

best_aic_ma <- Inf
best_order_ma <- c(0, 0)
best_model_ma <- NULL

for (q in q_values) {
  model <- tryCatch({
    arima(sales_ts, order = c(0, 0, q))  
  }, error = function(e) NULL)  
    
  if (!is.null(model)) {
    current_aic <- AIC(model)  
    
    if (current_aic < best_aic_ma) {
      best_aic_ma <- current_aic
      best_order_ma <- c(0, q)
      best_model_ma <- model
    }
  }
}

modelos <- c("ARMA", "ARIMA", "AR", "MA")
mejores_aic <- c(best_aic_arma, best_aic_arima, best_aic_ar, best_aic_ma)

comparacion_modelos <- data.frame(
  Modelo = modelos,
  AIC = mejores_aic
)

print("Comparación de modelos y sus AIC:")

## [1] "Comparación de modelos y sus AIC:"

print(comparacion_modelos)

##   Modelo      AIC
## 1   ARMA 1398.387
## 2  ARIMA 1354.534
## 3     AR 1401.696
## 4     MA 1402.149

mejor_modelo <- modelos[which.min(mejores_aic)]
mejor_aic <- min(mejores_aic)

print(paste("El mejor modelo es:", mejor_modelo))

## [1] "El mejor modelo es: ARIMA"

print(paste("El mejor AIC es:", mejor_aic))

## [1] "El mejor AIC es: 1354.53371359874"

if (mejor_modelo == "ARMA") {
  summary(best_model_arma)
} else if (mejor_modelo == "ARIMA") {
  summary(best_model_arima)
} else if (mejor_modelo == "AR") {
  summary(best_model_ar)
} else if (mejor_modelo == "MA") {
  summary(best_model_ma)
}

## 
## Call:
## arima(x = sales_ts, order = c(p, d, q))
## 
## Coefficients:
##           ar1      ar2      ma1      ma2      ma3     ma4     ma5
##       -0.7718  -0.5521  -0.5689  -0.3367  -0.8268  0.0062  0.7522
## s.e.   0.1690   0.1435   0.1924   0.2089   0.1744  0.1985  0.1731
## 
## sigma^2 estimated as 1.885e+11:  log likelihood = -669.27,  aic = 1352.53
## 
## Training set error measures:
##               ME RMSE MAE MPE MAPE
## Training set NaN  NaN NaN NaN  NaN

Forecast using ARIMA

modelo_arima <- Arima(sales_ts, order = c(2, 2, 5))


forecast_arima <- forecast(modelo_arima, h = 5)

autoplot(forecast_arima) +
  labs(title = "Pronóstico a 5 períodos con ARIMA(2,2,5)", 
       x = "Tiempo", 
       y = "Ventas (Cajas Unitarias)") +
  theme_minimal()

Time series model 2

Instructions: From the time series dataset, select 2 explanatory variables that might affect the sales unitboxes.

correlation_matrix <- cor(cocacola[, c("sales_unitboxes", 
                                       "consumer_sentiment", 
                                       "gdp_percapita", 
                                       "CPI", 
                                       "inflation_rate", 
                                       "unemp_rate", 
                                       "exchange_rate", 
                                       "max_temperature", 
                                       "Moving_Avg")], 
                          use = "complete.obs")


print(correlation_matrix)

##                    sales_unitboxes consumer_sentiment gdp_percapita         CPI
## sales_unitboxes         1.00000000         0.35524547    0.05786549  0.08807625
## consumer_sentiment      0.35524547         1.00000000   -0.52631149 -0.40446243
## gdp_percapita           0.05786549        -0.52631149    1.00000000  0.86618848
## CPI                     0.08807625        -0.40446243    0.86618848  1.00000000
## inflation_rate         -0.55212057        -0.52620078   -0.01727899  0.10332937
## unemp_rate              0.04231135         0.51952301   -0.76038173 -0.80824729
## exchange_rate          -0.06014764        -0.73031564    0.66028999  0.54067761
## max_temperature         0.68348780         0.03848755    0.37713494  0.17746354
## Moving_Avg              0.19129530        -0.20416077    0.73121504  0.65558384
##                    inflation_rate  unemp_rate exchange_rate max_temperature
## sales_unitboxes       -0.55212057  0.04231135   -0.06014764      0.68348780
## consumer_sentiment    -0.52620078  0.51952301   -0.73031564      0.03848755
## gdp_percapita         -0.01727899 -0.76038173    0.66028999      0.37713494
## CPI                    0.10332937 -0.80824729    0.54067761      0.17746354
## inflation_rate         1.00000000 -0.19810803    0.41019512     -0.55467320
## unemp_rate            -0.19810803  1.00000000   -0.63120125     -0.11731604
## exchange_rate          0.41019512 -0.63120125    1.00000000      0.13404524
## max_temperature       -0.55467320 -0.11731604    0.13404524      1.00000000
## Moving_Avg            -0.10970141 -0.51257546    0.52480755      0.25792091
##                    Moving_Avg
## sales_unitboxes     0.1912953
## consumer_sentiment -0.2041608
## gdp_percapita       0.7312150
## CPI                 0.6555838
## inflation_rate     -0.1097014
## unemp_rate         -0.5125755
## exchange_rate       0.5248076
## max_temperature     0.2579209
## Moving_Avg          1.0000000

sales_correlations <- correlation_matrix[1, -1]  

print(sales_correlations[order(abs(sales_correlations), decreasing = TRUE)])

##    max_temperature     inflation_rate consumer_sentiment         Moving_Avg 
##         0.68348780        -0.55212057         0.35524547         0.19129530 
##                CPI      exchange_rate      gdp_percapita         unemp_rate 
##         0.08807625        -0.06014764         0.05786549         0.04231135

max_temperature_diff <- diff(cocacola$max_temperature)
inflation_rate_diff <- diff(cocacola$inflation_rate)

Describe the hypothetical relationship / impact between each selected explanatory variable and the dependent variable “sales unit boxes”.

Temperature relationship with sales

It is extremely likely that our temperature variable has a positive correlation with our dependent variable. This is due to the fact that people are inclined towards buying Coca-Cola i warmer months, due to it’s refreshing nature and qualities. On the other hand, in colder months people are less icnlined to buy Coca-Cola, while products like coffee and tea might spike in popularity.

Inflation relationship with sales

There is likely a negative relationship between “inflation_rate” and sales unit boxes of Coca-Cola. When the inflation goes up, money loses value, and so people around the country lose buying power and effective income. For this reason, people are inclined to buy essential items, while dropping expenditures in things that can be avoided. One such example of this is Coca-Cola.

Include a time series plot that displays the selected variables’ performance over the time period.

aligned_data <- data.frame(
  tperiod = cocacola$tperiod[-1],  
  sales_unitboxes = cocacola$sales_unitboxes[-1],  
  max_temperature_diff = max_temperature_diff,  
  inflation_rate_diff = inflation_rate_diff  
)

ggplot() +
  geom_line(data = aligned_data, aes(x = tperiod, y = sales_unitboxes, color = "Ventas (Cajas Unitarias)"), size = 1) +
  
  geom_line(data = aligned_data, aes(x = tperiod, y = max_temperature_diff * 200000, color = "Temperatura Máxima (Diferenciada)"), size = 1, linetype = "dashed") +
  
  geom_line(data = aligned_data, aes(x = tperiod, y = inflation_rate_diff * 2000000, color = "Tasa de Inflación (Diferenciada)"), size = 1, linetype = "dotted") +
  
  scale_y_continuous(
    name = "Ventas (Cajas Unitarias)",
    sec.axis = sec_axis(~ . / 100000, name = "Temperatura (°C) y Tasa de Inflación Diferenciada (%)")
  ) +
  
  labs(title = "Ventas, Temperatura Máxima y Tasa de Inflación Diferenciada con Ejes Secundarios",
       x = "Tiempo", color = "") +
  
  theme_minimal() +
  theme(legend.position = "bottom")

Do the selected explanatory variables show serial autocorrelation and a stationary process?

acf_max_temperature_diff <- Acf(max_temperature_diff, main = "ACF de max_temperature (Diferenciada)")

acf_inflation_rate_diff <- Acf(inflation_rate_diff, main = "ACF de inflation_rate (Diferenciada)")

adf_test_max_temp_diff <- adf.test(max_temperature_diff, alternative = "stationary")
print(adf_test_max_temp_diff)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  max_temperature_diff
## Dickey-Fuller = -4.0323, Lag order = 3, p-value = 0.01639
## alternative hypothesis: stationary

adf_test_inflation_diff <- adf.test(inflation_rate_diff, alternative = "stationary")
print(adf_test_inflation_diff)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  inflation_rate_diff
## Dickey-Fuller = -3.895, Lag order = 3, p-value = 0.02207
## alternative hypothesis: stationary

Explanation of test results

Max Temperature: The p-value of lower than 0.05 tells us that our alternative hypothesis of the series being stationary is acceptable. The ADF also supports this statement, since it’s values are lower than critical (-3.43).

Inflation Rate: We can reject the null hypothesis, due to the fact that the p-value is lower than 0.05 (0.022). The ADF test statistic (-3.90) is also below critical values, allowing us to reject the null hypothesis of being not-stationary.

ADF Conclusion: The properties for both series are remaining constant over time, even after the application of differences. This makes them useful towards utilizing them in time series analysis technique involving stationarity,

Autocorrelation analysis after differentiation: In both cases we can see autocorrelation towards the start of the series (lag 1). This tells us that the data isn’t 100% free of autocorrelation, but this value decreases for all further lags. This tells us that the differencing had a positive impact on our model, since it heavily decreased the autocorrelation between the series after the first lag.

Briefly interpret the estimated VAR regression results. That is, is there a statistically significant relationship between the explanatory variable(s) and the dependent variable?

var_data <- data.frame(
  sales_unitboxes = cocacola$sales_unitboxes[-1],  
  inflation_rate_diff = inflation_rate_diff  
)

lag_selection <- VARselect(var_data, lag.max = 4, type = "const")
optimal_lag <- lag_selection$selection["AIC(n)"]  

VAR_model <- VAR(var_data, p = optimal_lag, type = "const")

summary(VAR_model)

## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: sales_unitboxes, inflation_rate_diff 
## Deterministic variables: const 
## Sample size: 44 
## Log Likelihood: -636.938 
## Roots of the characteristic polynomial:
## 0.7906 0.7906 0.7632 0.7632 0.7044 0.7044
## Call:
## VAR(y = var_data, p = optimal_lag, type = "const")
## 
## 
## Estimation results for equation sales_unitboxes: 
## ================================================ 
## sales_unitboxes = sales_unitboxes.l1 + inflation_rate_diff.l1 + sales_unitboxes.l2 + inflation_rate_diff.l2 + sales_unitboxes.l3 + inflation_rate_diff.l3 + const 
## 
##                          Estimate Std. Error t value Pr(>|t|)   
## sales_unitboxes.l1      4.718e-01  1.514e-01   3.117  0.00353 **
## inflation_rate_diff.l1 -5.059e+05  2.238e+05  -2.260  0.02978 * 
## sales_unitboxes.l2     -1.179e-02  1.726e-01  -0.068  0.94590   
## inflation_rate_diff.l2 -2.666e+05  2.408e+05  -1.107  0.27529   
## sales_unitboxes.l3      6.755e-02  1.649e-01   0.410  0.68445   
## inflation_rate_diff.l3 -5.693e+05  2.055e+05  -2.771  0.00870 **
## const                   3.105e+06  1.136e+06   2.733  0.00956 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 463100 on 37 degrees of freedom
## Multiple R-Squared: 0.4425,  Adjusted R-squared: 0.3521 
## F-statistic: 4.895 on 6 and 37 DF,  p-value: 0.0008983 
## 
## 
## Estimation results for equation inflation_rate_diff: 
## ==================================================== 
## inflation_rate_diff = sales_unitboxes.l1 + inflation_rate_diff.l1 + sales_unitboxes.l2 + inflation_rate_diff.l2 + sales_unitboxes.l3 + inflation_rate_diff.l3 + const 
## 
##                          Estimate Std. Error t value Pr(>|t|)    
## sales_unitboxes.l1      6.975e-08  9.946e-08   0.701  0.48752    
## inflation_rate_diff.l1 -6.279e-01  1.470e-01  -4.271  0.00013 ***
## sales_unitboxes.l2      1.302e-07  1.134e-07   1.148  0.25825    
## inflation_rate_diff.l2 -4.187e-01  1.582e-01  -2.647  0.01186 *  
## sales_unitboxes.l3      3.463e-07  1.083e-07   3.196  0.00285 ** 
## inflation_rate_diff.l3 -2.437e-01  1.350e-01  -1.805  0.07922 .  
## const                  -3.523e+00  7.463e-01  -4.721 3.33e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.3042 on 37 degrees of freedom
## Multiple R-Squared: 0.4472,  Adjusted R-squared: 0.3576 
## F-statistic: 4.989 on 6 and 37 DF,  p-value: 0.0007832 
## 
## 
## 
## Covariance matrix of residuals:
##                     sales_unitboxes inflation_rate_diff
## sales_unitboxes           2.145e+11          -4.109e+04
## inflation_rate_diff      -4.109e+04           9.256e-02
## 
## Correlation matrix of residuals:
##                     sales_unitboxes inflation_rate_diff
## sales_unitboxes              1.0000             -0.2917
## inflation_rate_diff         -0.2917              1.0000

forecast_steps <- 5
var_forecast <- predict(VAR_model, n.ahead = forecast_steps)

forecast_index <- seq(as.Date(tail(cocacola$tperiod, 1)), by = "month", length.out = forecast_steps)

forecast_df <- data.frame(
  Period = forecast_index,
  sales_unitboxes_pred = var_forecast$fcst$sales_unitboxes[, 1], 
  inflation_rate_pred = var_forecast$fcst$inflation_rate_diff[, 1]  
)

print(forecast_df)

##       Period sales_unitboxes_pred inflation_rate_pred
## 1 2018-12-01              6485870         -0.13676581
## 2 2019-01-01              6423873         -0.02262784
## 3 2019-02-01              6635569          0.14764207
## 4 2019-03-01              6607301         -0.02813486
## 5 2019-04-01              6565818         -0.01274115

There is a significant relationship between sales and both its own previous values (lag 1) and the inflation rate (lags 1 and 3). The inflation rate generally has a negative impact on sales, while past sales positively influence current sales.

The inflation rate is mainly affected by its own previous values, which have a strong negative impact. There is also some evidence that past sales influence inflation rate changes at lag 3.

Overall, the dependent variable is significantly influenced by its own past values and the inflation rate. Inflation tends to reduce sales, while past sales help to increase them. The inflation rate itself is influenced by its past values and past sales.

granger_test_1 <- causality(VAR_model, cause = "inflation_rate_diff")
print(granger_test_1$Granger)

## 
##  Granger causality H0: inflation_rate_diff do not Granger-cause
##  sales_unitboxes
## 
## data:  VAR object VAR_model
## F-Test = 3.7763, df1 = 3, df2 = 74, p-value = 0.01403

granger_test_2 <- causality(VAR_model, cause = "sales_unitboxes")
print(granger_test_2$Granger)

## 
##  Granger causality H0: sales_unitboxes do not Granger-cause
##  inflation_rate_diff
## 
## data:  VAR object VAR_model
## F-Test = 7.9594, df1 = 3, df2 = 74, p-value = 0.0001142

• Based on the selected VAR_Model, show the forecast of coca sales in the Guadalajara Metropolitan Area for the next 5 periods. Display the estimated forecast in a time series plot.

combined_data <- data.frame(
  Period = c(as.Date(cocacola$tperiod[-1]), forecast_df$Period), 
  sales_unitboxes = c(cocacola$sales_unitboxes[-1], rep(NA, forecast_steps)),  
  sales_unitboxes_pred = c(rep(NA, length(cocacola$sales_unitboxes[-1])), forecast_df$sales_unitboxes_pred)
)

ggplot() +
  geom_line(data = combined_data, aes(x = Period, y = sales_unitboxes, color = "Sales (Real)"), size = 1) +
  
  geom_line(data = combined_data, aes(x = Period, y = sales_unitboxes_pred, color = "Sales (Forecast)"), size = 1, linetype = "dashed") +
  
  geom_vline(xintercept = as.numeric(forecast_df$Period[1]), color = "yellow", linetype = "dashed", size = 1) +
  
  
  labs(title = "Monthly Sales and VAR forecast of Unit Boxes",
       x = "Date", y = "Sales (Unit Boxes)") +
  
 
  scale_color_manual(values = c("Sales (Real)" = "blue", "Sales (Forecast)" = "green")) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(legend.position = "bottom") +
  theme(legend.title = element_blank())

V. Conclusions and Recommendations

Conclusions

1. While sales_unitboxes Granger-cause inflation_rate_diff, the reverse is not true. This suggests that past values of sales_unitboxes provide significant information for predicting changes in inflation rate.

2. The rejection of the null hypothesis in the second Granger-causality test indicates that sales unit boxes have a predictive influence on inflation rates after differencing. This might suggest that changes in sales activity could be an indicator of inflationary trends.

3. Since the first Granger-causality test fails to reject the null hypothesis , inflation (after differencing) does not appear to Granger-cause sales unit boxes. This could imply that fluctuations in inflation may not be directly influencing sales activity

4. The covariance matrix of residuals shows a weak negative correlation (-0.2917) between sales_unitboxes and inflation_rate_diff. This suggests that, while the residuals of both variables are interconnected, they do not exhibit a strong linear relationship, implying that other factors may influence the dynamics between sales and inflation beyond what is captured by their past values.

5. The VAR forecast begins in late 2018, just after a moderate upward trend in real sales. The forecast shows a relatively stable prediction for future sales, with values slightly above the final real sales point. This indicates that the model expects a stabilization of sales based on recent patterns, projecting moderate growth or consistency in the near future.

6. The forecasted section of the graph maintains the clear seasonal patterns seen in historical data, which suggests that these patterns will persist. Coca-Cola can use this information to prepare for peak demand periods, ensuring inventory and marketing efforts are aligned with these predictable sales surges.

Recommendations

1. Adjust prices for Coca-Cola products based on the weather. For example, offer discounts or promotions on hotter days to encourage more sales. In colder months, you might focus on bundling with other products.

2. Implement a pricing strategy that adjusts to inflationary pressures. For example, offering the product in smaller containers, might make it easier for customers to buy it, mainly in countries like Mexico, where the income is low in many families.

3. We could also consider using machine learning models like neural networks, which can handle time series data and complex non-linear relationships. By feeding these models with relevant time series features, we can improve the accuracy of the forecasts we made.