Granger Causality

# ============================================================================
# GRANGER CAUSALITY TEST: Does Advertising Granger-Cause Sales?
# ============================================================================
# Question: Does advertising expenditure (X) Granger-cause sales revenue (Y)?
# Data: 156 weeks (3 years) of weekly advertising and sales data
# Method: Granger causality test via F-test on VAR model
# ============================================================================
# CRITICAL: Granger causality ≠ true causation
#           It tests if X's past helps PREDICT Y, not if X CAUSES Y
# ============================================================================

# Load required packages
library(lmtest)        # grangertest()
library(vars)          # VAR(), causality()
library(tseries)       # adf.test()
library(ggplot2)       # visualization
library(gridExtra)     # multiple plots
library(forecast)      # auto.arima for lag selection

set.seed(123)

# ============================================================================
# 1. DATA GENERATION: Advertising and Sales with Known Causal Relationship
# ============================================================================

cat("============================================================================\n")
cat("GRANGER CAUSALITY TEST: Advertising → Sales\n")
cat("============================================================================\n\n")

# Generate 156 weeks of data (3 years)
n <- 156
time_index <- 1:n

# Generate advertising expenditure (X) - AR(1) process with trend
adv_innovation <- rnorm(n, 0, 5)
advertising <- numeric(n)
advertising[1] <- 50
for (t in 2:n) {
  advertising[t] <- 40 + 0.1*t + 0.5*advertising[t-1] + adv_innovation[t]
}

# Generate sales (Y) - depends on LAGGED advertising (2-week lag) + own dynamics
sales_innovation <- rnorm(n, 0, 10)
sales <- numeric(n)
sales[1] <- 100
sales[2] <- 100
for (t in 3:n) {
  # Sales depends on: own lag, advertising from 2 weeks ago, and noise
  sales[t] <- 80 + 0.3*sales[t-1] + 0.8*advertising[t-2] + sales_innovation[t]
}

# Create data frame
data <- data.frame(
  week = time_index,
  advertising = advertising,
  sales = sales
)

cat("Data generated: n =", n, "weeks\n")
cat("TRUE relationship: Sales[t] = f(Sales[t-1], Advertising[t-2]) + noise\n")
cat("  → Advertising DOES Granger-cause Sales(2-week lag)\n\n")

cat("Summary statistics:\n")
cat("Advertising: Mean=", round(mean(data$advertising), 2), 
    ", SD=", round(sd(data$advertising), 2), "\n")
cat("Sales: Mean=", round(mean(data$sales), 2), 
    ", SD=", round(sd(data$sales), 2), "\n\n")

# ============================================================================
# 2. EXPLORATORY VISUALIZATION
# ============================================================================

cat("=== EXPLORATORY ANALYSIS ===\n\n")

# Time series plots
p1 <- ggplot(data, aes(x=week, y=advertising)) +
  geom_line(color="steelblue", size=1) +
  labs(title="Advertising Expenditure over Time", x="Week", y="Advertising($1000s)") +
  theme_minimal()

p2 <- ggplot(data, aes(x=week, y=sales)) +
  geom_line(color="darkgreen", size=1) +
  labs(title="Sales Revenue over Time", x="Week", y="Sales($1000s)") +
  theme_minimal()

# Scatter plot: Sales vs Lagged Advertising
data$adv_lag2 <- c(NA, NA, data$advertising[1:(n-2)])

p3 <- ggplot(data, aes(x=adv_lag2, y=sales)) +
  geom_point(alpha=0.5, color="darkred") +
  geom_smooth(method="lm", se=TRUE, color="blue") +
  labs(title="Sales vs Advertising(2-week lag)", 
       x="Advertising[t-2]", y="Sales[t]") +
  theme_minimal()

# Cross-correlation function (CCF)
par(mfrow=c(1,1))
ccf_result <- ccf(data$advertising, data$sales, lag.max=20, 
                  main="Cross-Correlation: Advertising vs Sales",
                  ylab="CCF", xlab="Lag(negative=advertising leads sales)")

cat("Cross-correlation function (CCF) computed\n")
cat("  Negative lags: advertising LEADS sales\n")
cat("  Positive lags: sales LEADS advertising\n")
cat("  Peak at negative lag suggests advertising predicts future sales\n\n")

# Find peak CCF
peak_lag <- which.max(abs(ccf_result$acf)) - 21  # Adjust for 0-lag at position 21
cat("Peak CCF at lag:", peak_lag, "\n")
if (peak_lag < 0) {
  cat("  → Advertising leads sales by", abs(peak_lag), "weeks(consistent with Granger causality)\n\n")
}

# ============================================================================
# 3. STATIONARITY TESTING
# ============================================================================

cat("=== STATIONARITY TESTS ===\n\n")

# Function to test stationarity
test_stationarity <- function(series, name) {
  cat(name, ":\n")
  
  # ADF test
  adf <- adf.test(series)
  cat("  ADF test: statistic=", round(adf$statistic, 4), 
      ", p-value=", round(adf$p.value, 4), "\n")
  if (adf$p.value < 0.05) {
    cat("    → STATIONARY(p < 0.05)\n")
  } else {
    cat("    → NON-STATIONARY(p ≥ 0.05)\n")
  }
  
  # KPSS test
  kpss <- kpss.test(series, null="Trend")
  cat("  KPSS test: statistic=", round(kpss$statistic, 4), 
      ", p-value=", round(kpss$p.value, 4), "\n")
  if (kpss$p.value >= 0.05) {
    cat("    → STATIONARY(p ≥ 0.05)\n")
  } else {
    cat("    → NON-STATIONARY(p < 0.05)\n")
  }
  cat("\n")
  
  return(adf$p.value < 0.05 & kpss$p.value >= 0.05)
}

# Test both series
adv_stationary <- test_stationarity(data$advertising, "Advertising")
sales_stationary <- test_stationarity(data$sales, "Sales")

if (!adv_stationary | !sales_stationary) {
  cat("WARNING: At least one series non-stationary\n")
  cat("Action: Apply first differencing\n\n")
  
  # Difference both series
  data$adv_diff <- c(NA, diff(data$advertising))
  data$sales_diff <- c(NA, diff(data$sales))
  
  cat("Testing differenced series:\n")
  test_stationarity(na.omit(data$adv_diff), "Advertising(differenced)")
  test_stationarity(na.omit(data$sales_diff), "Sales(differenced)")
  
  # Use differenced series
  use_diff <- TRUE
  cat("Proceeding with DIFFERENCED series for Granger test\n")
  cat("NOTE: Interpretation changes to 'changes in X cause changes in Y'\n\n")
} else {
  use_diff <- FALSE
  cat("Both series stationary → proceed with LEVELS\n\n")
}

# ============================================================================
# 4. LAG ORDER SELECTION
# ============================================================================

cat("=== LAG ORDER SELECTION ===\n\n")

if (use_diff) {
  # VAR on differenced data (remove NA from differencing)
  var_data <- data.frame(
    advertising = data$adv_diff,
    sales = data$sales_diff
  )
  var_data <- na.omit(var_data)
} else {
  var_data <- data.frame(
    advertising = data$advertising,
    sales = data$sales
  )
}

# Select lag order using information criteria
lag_select <- VARselect(var_data, lag.max=12, type="const")

cat("Information criteria for VAR lag selection:\n")
print(lag_select$selection)
cat("\n")

aic_lag <- lag_select$selection["AIC(n)"]
bic_lag <- lag_select$selection["SC(n)"]  # SC = Schwarz = BIC

cat("Selected lags:\n")
cat("  AIC recommends p =", aic_lag, "\n")
cat("  BIC recommends p =", bic_lag, "(more conservative)\n\n")

# Use BIC lag (more conservative, avoids overfitting)
optimal_lag <- bic_lag
cat("Using p =", optimal_lag, "lags(BIC criterion)\n\n")

# ============================================================================
# 5. GRANGER CAUSALITY TEST: Does Advertising → Sales?
# ============================================================================

cat("=== GRANGER CAUSALITY TEST: Advertising → Sales ===\n\n")

# Method 1: Using grangertest() from lmtest
cat("Method 1: grangertest() from lmtest package\n")
cat("---------------------------------------------\n")

granger_adv_sales <- grangertest(sales ~ advertising, 
                                  order=optimal_lag, data=var_data)
print(granger_adv_sales)
cat("\n")

if (granger_adv_sales$`Pr(>F)`[2] < 0.05) {
  cat("RESULT: REJECT H₀ (p < 0.05)\n")
  cat("  → Advertising DOES Granger-cause Sales\n")
  cat("  → Past advertising helps PREDICT future sales\n")
} else {
  cat("RESULT: FAIL TO REJECT H₀ (p ≥ 0.05)\n")
  cat("  → Advertising does NOT Granger-cause Sales\n")
  cat("  → Past advertising does not help predict sales\n")
}
cat("\n")

# Method 2: Using VAR() and causality() from vars package
cat("Method 2: VAR model with causality() test\n")
cat("------------------------------------------\n")

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

# Granger causality test via vars::causality()
causality_test <- causality(var_model, cause="advertising")

cat("\nGranger causality(F-test): Advertising → Sales\n")
print(causality_test$Granger)

cat("\nInterpretation:\n")
if (causality_test$Granger$p.value < 0.05) {
  cat("  p-value =", round(causality_test$Granger$p.value, 4), "< 0.05\n")
  cat("  → SIGNIFICANT: Advertising Granger-causes Sales\n\n")
} else {
  cat("  p-value =", round(causality_test$Granger$p.value, 4), "≥ 0.05\n")
  cat("  → NOT SIGNIFICANT: No Granger causality detected\n\n")
}

# ============================================================================
# 6. BI-DIRECTIONAL TESTING: Does Sales → Advertising?
# ============================================================================

cat("=== BI-DIRECTIONAL TEST: Sales → Advertising ===\n\n")

granger_sales_adv <- grangertest(advertising ~ sales, 
                                  order=optimal_lag, data=var_data)
print(granger_sales_adv)
cat("\n")

if (granger_sales_adv$`Pr(>F)`[2] < 0.05) {
  cat("RESULT: Sales DOES Granger-cause Advertising\n")
  cat("  → Past sales help predict future advertising\n")
  cat("  → BIDIRECTIONAL causality(feedback loop)\n")
} else {
  cat("RESULT: Sales does NOT Granger-cause Advertising\n")
  cat("  → UNIDIRECTIONAL causality(Advertising → Sales only)\n")
}
cat("\n")

# Summary of causality direction
cat("CAUSALITY SUMMARY:\n")
cat("------------------\n")
adv_to_sales <- granger_adv_sales$`Pr(>F)`[2] < 0.05
sales_to_adv <- granger_sales_adv$`Pr(>F)`[2] < 0.05

if (adv_to_sales & sales_to_adv) {
  cat("  Direction: BIDIRECTIONAL(X ↔ Y)\n")
  cat("  Advertising ↔ Sales: Feedback loop\n")
} else if (adv_to_sales & !sales_to_adv) {
  cat("  Direction: UNIDIRECTIONAL(X → Y)\n")
  cat("  Advertising → Sales only\n")
} else if (!adv_to_sales & sales_to_adv) {
  cat("  Direction: REVERSE(Y → X)\n")
  cat("  Sales → Advertising only\n")
} else {
  cat("  Direction: NONE\n")
  cat("  No Granger causality detected in either direction\n")
}
cat("\n")

# ============================================================================
# 7. VAR MODEL COEFFICIENTS - Show Which Lags Matter
# ============================================================================

cat("=== VAR MODEL COEFFICIENTS ===\n\n")

cat("Sales equation(showing which lags of Advertising predict Sales):\n")
cat("---------------------------------------------------------------\n")
summary(var_model)$varresult$sales

cat("\nKey coefficients to examine:\n")
sales_coefs <- coef(var_model$varresult$sales)
adv_coefs <- sales_coefs[grep("advertising", names(sales_coefs))]

cat("\nLagged Advertising coefficients in Sales equation:\n")
for (i in 1:length(adv_coefs)) {
  coef_val <- adv_coefs[i]
  coef_name <- names(adv_coefs)[i]
  sig <- abs(coef_val / sqrt(vcov(var_model$varresult$sales)[i,i])) > 1.96
  
  cat(" ", coef_name, ": ", round(coef_val, 4))
  if (sig) cat(" ***")
  cat("\n")
}
cat("*** = significant at 5% level(|t| > 1.96)\n\n")

# ============================================================================
# 8. IMPULSE RESPONSE FUNCTIONS (IRF)
# ============================================================================

cat("=== IMPULSE RESPONSE FUNCTIONS ===\n\n")

# IRF: How does Sales respond to shock in Advertising?
irf_result <- irf(var_model, impulse="advertising", response="sales", 
                   n.ahead=20, boot=TRUE, runs=1000)

cat("Computing IRF: Effect of Advertising shock on Sales over 20 periods\n")
cat("Using bootstrap(1000 runs) for confidence intervals\n\n")

plot(irf_result, main="Impulse Response: Sales response to Advertising shock")

cat("\nInterpretation:\n")
cat("  - Shows dynamic effect of one-time advertising shock on sales\n")
cat("  - Y-axis: change in sales\n")
cat("  - X-axis: periods after shock\n")
cat("  - If response is positive and significant: advertising shock increases sales\n\n")

# ============================================================================
# 9. FORECAST ERROR VARIANCE DECOMPOSITION (FEVD)
# ============================================================================

cat("=== FORECAST ERROR VARIANCE DECOMPOSITION ===\n\n")

# FEVD: What % of Sales forecast error is due to Advertising shocks?
fevd_result <- fevd(var_model, n.ahead=20)

cat("FEVD for Sales(% of forecast error variance explained by each variable):\n")
cat("-------------------------------------------------------------------------\n")
print(round(fevd_result$sales * 100, 2))
cat("\n")

cat("Interpretation:\n")
adv_contribution <- fevd_result$sales[10, "advertising"]  # 10 periods ahead
cat("  At 10-period horizon, Advertising explains", 
    round(adv_contribution*100, 1), "% of Sales forecast error variance\n")
if (adv_contribution > 0.1) {
  cat("  → Advertising is IMPORTANT for Sales prediction\n")
} else {
  cat("  → Advertising has WEAK explanatory power(other factors dominate)\n")
}
cat("\n")

# ============================================================================
# 10. RESIDUAL DIAGNOSTICS
# ============================================================================

cat("=== RESIDUAL DIAGNOSTICS ===\n\n")

# Portmanteau test for residual autocorrelation
serial_test <- serial.test(var_model, lags.pt=16, type="PT.asymptotic")

cat("Portmanteau test for residual autocorrelation:\n")
print(serial_test)
cat("\n")

if (serial_test$serial$p.value > 0.05) {
  cat("GOOD: Residuals are white noise(no autocorrelation)\n")
  cat("  → VAR model adequately captures dynamics\n")
} else {
  cat("WARNING: Residuals show autocorrelation\n")
  cat("  → May need to increase lag order\n")
}
cat("\n")

# Plot residuals
par(mfrow=c(2,2))
plot(var_model$varresult$sales$residuals, type="l", 
     main="Sales Equation Residuals", ylab="Residuals")
abline(h=0, col="red", lty=2)

acf(var_model$varresult$sales$residuals, main="ACF: Sales Residuals")

plot(var_model$varresult$advertising$residuals, type="l", 
     main="Advertising Equation Residuals", ylab="Residuals")
abline(h=0, col="red", lty=2)

acf(var_model$varresult$advertising$residuals, main="ACF: Advertising Residuals")
par(mfrow=c(1,1))

# ============================================================================
# 11. ROBUSTNESS: Test at Multiple Lags
# ============================================================================

cat("=== ROBUSTNESS CHECK: Multiple Lag Orders ===\n\n")

lag_range <- 1:8
p_values <- numeric(length(lag_range))

for (i in seq_along(lag_range)) {
  p <- lag_range[i]
  gt <- grangertest(sales ~ advertising, order=p, data=var_data)
  p_values[i] <- gt$`Pr(>F)`[2]
}

robustness_df <- data.frame(
  Lag = lag_range,
  P_value = round(p_values, 4),
  Significant = ifelse(p_values < 0.05, "Yes", "No")
)

cat("Granger test p-values across different lag orders:\n")
print(robustness_df, row.names=FALSE)
cat("\n")

# Plot p-values across lags
par(mfrow=c(1,1))
plot(lag_range, p_values, type="b", pch=19, col="steelblue",
     xlab="Lag order(p)", ylab="p-value",
     main="Granger Causality p-values across Lag Orders",
     ylim=c(0, max(p_values, 0.1)))
abline(h=0.05, col="red", lty=2, lwd=2)
text(max(lag_range)*0.7, 0.06, "α = 0.05", col="red")

cat("Interpretation:\n")
sig_lags <- sum(p_values < 0.05)
cat(" ", sig_lags, "out of", length(lag_range), "lag orders show significant causality\n")
if (sig_lags >= length(lag_range)/2) {
  cat("  → ROBUST: Granger causality holds across multiple lag specifications\n")
} else {
  cat("  → WEAK: Causality sensitive to lag choice(may be spurious)\n")
}
cat("\n")

# ============================================================================
# 12. CRITICAL INTERPRETATION & WARNINGS
# ============================================================================

cat("============================================================================\n")
cat("CRITICAL INTERPRETATION\n")
cat("============================================================================\n\n")

cat("WHAT WE FOUND:\n")
if (granger_adv_sales$`Pr(>F)`[2] < 0.05) {
  cat("  ✓ Advertising DOES Granger-cause Sales(p < 0.05)\n")
  cat("  ✓ Past advertising expenditure helps PREDICT future sales\n")
  cat("  ✓ Including advertising lags improves sales forecasts\n\n")
} else {
  cat("  ✗ No Granger causality detected(p ≥ 0.05)\n")
  cat("  ✗ Past advertising does not improve sales predictions\n\n")
}

cat("WHAT THIS DOES NOT MEAN:\n")
cat("  ✗ Does NOT prove advertising CAUSES sales\n")
cat("  ✗ Does NOT rule out reverse causation(sales → advertising)\n")
cat("  ✗ Does NOT account for confounders(economy, competition, etc.)\n")
cat("  ✗ Does NOT imply increasing advertising will increase sales\n\n")

cat("ALTERNATIVE EXPLANATIONS:\n")
cat("  1. Reverse causation: High sales → increased advertising budget\n")
cat("  2. Common cause: Economic growth drives both advertising and sales\n")
cat("  3. Spurious correlation: Non-stationarity creates false relationship\n")
cat("  4. Omitted variables: Competitor actions, seasonality, etc.\n\n")

cat("GRANGER CAUSALITY = PREDICTIVE PRECEDENCE, NOT CAUSATION\n\n")

cat("============================================================================\n")
cat("GRANGER CAUSALITY TEST COMPLETE\n")
cat("============================================================================\n")

# ============================================================================
# END OF GRANGER CAUSALITY EXAMPLE
# ============================================================================