1. Load Packages & Data

Make sure your CSV file is in the same folder as this .Rmd file and update the filename in read_csv() if needed.

library(tidyverse)
library(forecast)
library(tseries)
library(lmtest)
library(broom)

# Load dataset (update the path/file name if needed)
data_path <- "C:/Users/zxu3/Documents/R/textanalysis/Avocados/Nov_2026/data3_joint.csv"

df <- read_csv(data_path, show_col_types = FALSE) %>%
  mutate(date = as.Date(as.character(date), format = "%m/%d/%Y")) %>%
  arrange(time_idx)

# Check dimensions and missing weeks
glimpse(df)
## Rows: 256
## Columns: 42
## $ week                <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,…
## $ date                <date> 2015-01-04, 2015-01-11, 2015-01-18, 2015-01-25, 2…
## $ total_value         <dbl> 14821625, 14628353, 14907504, 14378714, 19756206, …
## $ id                  <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,…
## $ sentiment           <dbl> 6.666667, 16.500000, 10.000000, 6.500000, 11.00000…
## $ trend_value         <dbl> 48, 50, 48, 51, 56, 42, 39, 40, 43, 44, 44, 44, 47…
## $ log_total_value     <dbl> 16.51160, 16.49847, 16.51738, 16.48126, 16.79898, …
## $ mean_sentiment      <dbl> 6.666667, 16.500000, 10.000000, 6.500000, 11.00000…
## $ n_tweets            <dbl> 489, 144, 174, 231, 130, 189, 87, 62, 66, 81, 69, …
## $ total_volume        <dbl> 15540212, 14274341, 13970970, 13147070, 21473717, …
## $ total_revenue_check <dbl> 14821625, 14628353, 14907504, 14378714, 19756206, …
## $ avg_price_weighted  <dbl> 0.9537595, 1.0248006, 1.0670342, 1.0936820, 0.9200…
## $ organic_volume      <dbl> 298492.6, 334774.2, 366101.2, 264471.5, 369524.0, …
## $ conventional_volume <dbl> 15241720, 13939566, 13604869, 12882599, 21104193, …
## $ organic_share       <dbl> 0.01920776, 0.02345287, 0.02620442, 0.02011639, 0.…
## $ log_total_volume    <dbl> 16.55894, 16.47397, 16.45249, 16.39171, 16.88234, …
## $ log_avg_price       <dbl> -0.04734375, 0.02449808, 0.06488305, 0.08955001, -…
## $ super_bowl_week     <dbl> 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ cinco_de_mayo_week  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ memorial_day_week   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ july4_week          <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ labor_day_week      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ summer_peak         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ holiday_season      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ time_idx            <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,…
## $ sin1                <dbl> 0.06875541, 0.18809942, 0.30471923, 0.41692588, 0.…
## $ cos1                <dbl> 0.9976335, 0.9821500, 0.9524422, 0.9089405, 0.8522…
## $ sin2                <dbl> 0.1371854, 0.3694837, 0.5804549, 0.7579216, 0.8916…
## $ cos2                <dbl> 0.990545388, 0.929237218, 0.814292378, 0.652345619…
## $ post                <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ time_c              <dbl> -138, -137, -136, -135, -134, -133, -132, -131, -1…
## $ post_time           <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ evt_1_20161019      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ evt_2_20170518      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ evt_3_20170827      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ evt_10_20171209     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ evt_4_20180305      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ evt_5_20180901      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ evt_6_20180930      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ evt_7_20181130      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ evt_8_20190530      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ evt_9_20190607      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…

2. Define Dependent Variable & External Regressors (xreg)

We use log_total_value as the target variable (so coefficients on the 10 event variables can be interpreted as approximate percentage changes). You can easily swap this to log_total_volume or log_avg_price.

Your dataset already contains: 1. 10 Interruption Dummies: evt_1_20161019 to evt_10_20171209 2. Seasonality (Fourier Terms): sin1, cos1, sin2, cos2 (capturing the 52-week annual wave) 3. Holiday Spikes: super_bowl_week, cinco_de_mayo_week, memorial_day_week, july4_week, labor_day_week, holiday_season

# Target series (change to log_total_volume or log_avg_price if desired)
y <- ts(df$log_total_value, frequency = 52)

# Select the 10 interruptions + seasonality + holiday controls
xreg_vars <- df %>%
  select(
    # 10 Interruptions (ordered chronologically)
    evt_1_20161019,
    evt_2_20170518,
    evt_3_20170827,
    evt_10_20171209,
    evt_4_20180305,
    evt_5_20180901,
    evt_6_20180930,
    evt_7_20181130,
    evt_8_20190530,
    evt_9_20190607,
    
    # Annual Fourier Seasonality
    sin1, cos1, sin2, cos2,
    
    # Holiday Controls
    super_bowl_week,
    cinco_de_mayo_week,
    memorial_day_week,
    july4_week,
    labor_day_week,
    holiday_season
  ) %>%
  as.matrix()

3. Exploratory Plot & Stationarity Check

# Plot time series with the 10 interruption dates marked
event_dates <- as.Date(c(
  "2016-10-19", "2017-05-18", "2017-08-27", "2017-12-09", "2018-03-05",
  "2018-09-01", "2018-09-30", "2018-11-30", "2019-05-30", "2019-06-07"
))

ggplot(df, aes(x = date, y = log_total_value)) +
  geom_line(color = "#2c3e50", linewidth = 0.7) +
  geom_vline(xintercept = event_dates, linetype = "dashed", color = "#e74c3c", alpha = 0.7) +
  labs(
    title = "Weekly Log Total Value with 10 Interruption Points",
    subtitle = "Dashed red lines indicate the start of each event dummy",
    x = "Date",
    y = "Log Total Value"
  ) +
  theme_minimal()

# Augmented Dickey-Fuller test on raw series
adf.test(df$log_total_value)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  df$log_total_value
## Dickey-Fuller = -2.849, Lag order = 6, p-value = 0.2183
## alternative hypothesis: stationary

4. Fit ARIMAX Model (Automatic Lag Selection)

Because the Fourier terms (sin1–cos2) and holiday flags already absorb the 52-week seasonality, we set seasonal = FALSE inside auto.arima(). This prevents overfitting and lets auto.arima() select the best non-seasonal lags (p, d, q) via AICc.

fit <- auto.arima(
  y,
  xreg = xreg_vars,
  seasonal = FALSE,
  stepwise = FALSE,
  approximation = FALSE
)

summary(fit)
## Series: y 
## Regression with ARIMA(5,0,0) errors 
## 
## Coefficients:
##          ar1     ar2     ar3     ar4     ar5  intercept  evt_1_20161019
##       0.5272  0.0295  0.0526  0.1401  0.2128    16.7156         -0.0309
## s.e.  0.0664  0.0724  0.0705  0.0706  0.0660     0.0776          0.0443
##       evt_2_20170518  evt_3_20170827  evt_10_20171209  evt_4_20180305
##               0.0072          0.0402           0.0334          0.0094
## s.e.          0.0415          0.0408           0.0399          0.0392
##       evt_5_20180901  evt_6_20180930  evt_7_20181130  evt_8_20190530
##               0.0422         -0.0200          0.1151          0.0468
## s.e.          0.0403          0.0401          0.0413          0.0484
##       evt_9_20190607    sin1     cos1   sin2    cos2  super_bowl_week
##              -0.0258  0.0270  -0.0888  6e-04  0.0046           0.2394
## s.e.          0.0472  0.0155   0.0166  9e-03  0.0087           0.0189
##       cinco_de_mayo_week  memorial_day_week  july4_week  labor_day_week
##                   0.1429             0.0481      0.0725          0.0324
## s.e.              0.0191             0.0193      0.0189          0.0192
##       holiday_season
##              -0.1110
## s.e.          0.0167
## 
## sigma^2 = 0.00254:  log likelihood = 414.17
## AIC=-774.35   AICc=-767.71   BIC=-678.63
## 
## Training set error measures:
##                       ME       RMSE        MAE        MPE      MAPE      MASE
## Training set 0.002098209 0.04776657 0.03415661 0.01173555 0.2040929 0.2865559
##                     ACF1
## Training set -0.02005389

5. Statistical Significance & Impact of the 10 Interruptions

Since the dependent variable is logged (log_total_value), the exact percentage impact of each step interruption is calculated as \((\exp(\beta) - 1) \times 100\).

# Extract coefficients, standard errors, z-values, and p-values
coef_table <- coeftest(fit) %>%
  tidy() %>%
  mutate(
    pct_change = ifelse(grepl("^evt_", term), (exp(estimate) - 1) * 100, NA),
    significant = ifelse(p.value < 0.05, "Yes (*)", "No")
  )

# Display only the 10 interruption results cleanly
coef_table %>%
  filter(grepl("^evt_", term)) %>%
  knitr::kable(
    digits = 4,
    col.names = c("Interruption Variable", "Estimate (Log)", "Std. Error", "z-value", "p-value", "Est. % Change", "Significant (p < 0.05)")
  )
Interruption Variable Estimate (Log) Std. Error z-value p-value Est. % Change Significant (p < 0.05)
evt_1_20161019 -0.0309 0.0443 -0.6961 0.4864 -3.0383 No
evt_2_20170518 0.0072 0.0415 0.1743 0.8616 0.7262 No
evt_3_20170827 0.0402 0.0408 0.9845 0.3248 4.1026 No
evt_10_20171209 0.0334 0.0399 0.8361 0.4031 3.3938 No
evt_4_20180305 0.0094 0.0392 0.2406 0.8099 0.9464 No
evt_5_20180901 0.0422 0.0403 1.0474 0.2949 4.3108 No
evt_6_20180930 -0.0200 0.0401 -0.4988 0.6179 -1.9782 No
evt_7_20181130 0.1151 0.0413 2.7846 0.0054 12.2000 Yes (*)
evt_8_20190530 0.0468 0.0484 0.9674 0.3333 4.7896 No
evt_9_20190607 -0.0258 0.0472 -0.5468 0.5845 -2.5489 No

6. Residual Diagnostics (Checking Lag Adequacy)

If the selected (p, d, q) lags properly capture the time-series autocorrelation: * The ACF plot should show no significant spikes outside the blue dashed lines. * The Ljung-Box test p-value should be > 0.05 (confirming residuals behave like white noise).

checkresiduals(fit)

## 
##  Ljung-Box test
## 
## data:  Residuals from Regression with ARIMA(5,0,0) errors
## Q* = 73.719, df = 46, p-value = 0.00585
## 
## Model df: 5.   Total lags used: 51

7. Actual vs. Fitted Values

df_results <- df %>%
  mutate(
    fitted = as.numeric(fitted(fit)),
    residuals = as.numeric(residuals(fit))
  )

ggplot(df_results, aes(x = date)) +
  geom_line(aes(y = log_total_value, color = "Actual"), linewidth = 0.7) +
  geom_line(aes(y = fitted, color = "Fitted ARIMAX"), linewidth = 0.7, linetype = "solid", alpha = 0.85) +
  scale_color_manual(values = c("Actual" = "#2c3e50", "Fitted ARIMAX" = "#27ae60")) +
  labs(
    title = "Actual vs. Fitted Log Total Value",
    x = "Date",
    y = "Log Total Value",
    color = "Series"
  ) +
  theme_minimal()