基于R的ARIMAX+Fourier双销售序列关联影响分析及VAR方法咨询

Hey there! Let's work through your problem step by step—you're trying to understand how Series A would react if Series B hits its annual sales target, and you've been stuck finding solid VAR examples in R. First, let's cover why VAR makes sense here, then walk through a hands-on code example, and address the constraints you're dealing with (limited 2-year weekly data, seasonal trends, differing holiday responses).

Why VAR for Your Use Case? #

Unlike ARIMAX (which focuses on one-way effects from predictors to a single outcome), VAR models are built to capture dynamic, bidirectional relationships between correlated time series. This is perfect for your goal: you can simulate how a specific path for Series B (one that hits its annual target) would ripple through to Series A, while accounting for their shared trends, seasonality, and differing holiday behaviors.

Key Pre-Work for Your Data #

Two years of weekly data gives you ~104 observations—enough for a simple VAR, but we need to avoid overcomplicating the model. Here's how to handle your constraints:

• Seasonality: Stick with the Fourier terms you're already using, but feed them as exogenous variables into the VAR (instead of pre-processing the series). This preserves trend information while accounting for seasonality without overfitting.
• Holiday Differences: Create separate dummy variables for each key holiday (e.g., christmas_dummy, thanksgiving_dummy) instead of a single "holiday" flag. The VAR will automatically estimate how each series responds differently to each holiday.
• Stationarity: VAR works best with stationary series. Use the ADF test (tseries::adf.test()) to check—if your series are non-stationary, you can either difference them or include a linear trend term in the VAR (better for preserving trend-related context for annual targets).

Hands-On R VAR Example #

Let's turn this into actionable code, using simulated data that matches your scenario (adjust with your actual data):

# Load required packages
library(vars)
library(forecast)
library(dplyr)
library(lubridate)

# ----------------------
# Step 1: Prepare your data (replace with your actual dataset)
# ----------------------
set.seed(123) # For reproducibility
date_seq <- seq.Date(as.Date("2022-01-01"), as.Date("2023-12-31"), by = "week")

# Simulate Series A & B with shared trend + seasonality, plus noise
trend <- seq_along(date_seq) * 0.4
seasonal <- 8 * sin(2 * pi * seq_along(date_seq)/52) # Weekly seasonality (52 weeks/year)
A <- trend + seasonal + rnorm(length(date_seq), 0, 1.5)
B <- trend + seasonal * 1.3 + rnorm(length(date_seq), 0, 1.8) # B has stronger seasonality

# Add holiday dummies (customize to your actual holidays)
holiday_dummies <- data.frame(
  christmas = ifelse(week(date_seq) %in% c(51,52), 1, 0),
  thanksgiving = ifelse(week(date_seq) == 45, 1, 0)
)

# Combine into a dataframe
df <- bind_cols(
  data.frame(date = date_seq, A = A, B = B),
  holiday_dummies
)

# Convert to time series object (critical for VAR)
ts_data <- ts(df[, c("A", "B")], frequency = 52, start = c(2022, 1))

# ----------------------
# Step 2: Add Fourier terms for seasonality (exogenous variable)
# ----------------------
# Use K=2 for 2-year data (avoids overfitting—test K=1 or 3 if needed)
fourier_terms <- fourier(ts_data, K = 2)

# Combine all exogenous variables (Fourier + holidays)
exog_vars <- cbind(fourier_terms, df[, c("christmas", "thanksgiving")])
colnames(exog_vars) <- c("S1", "C1", "S2", "C2", "christmas", "thanksgiving")

# ----------------------
# Step 3: Select optimal lag order (critical for small samples)
# ----------------------
# Test up to 8 lags (weekly data—lags represent weeks)
var_lag_select <- VARselect(ts_data, lag.max = 8, exogen = exog_vars)
best_lag <- var_lag_select$selection["BIC(n)"] # BIC is more conservative for small data

# ----------------------
# Step 4: Fit the VAR model
# ----------------------
var_model <- VAR(
  y = ts_data,
  p = best_lag,
  exogen = exog_vars,
  type = "trend" # Include linear trend since your series have clear trends
)

# Check model diagnostics (ensure residuals are white noise)
serial.test(var_model) # If p-value > 0.05, residuals are uncorrelated (good!)

# ----------------------
# Step 5: Simulate "Series B hits annual target" scenario
# ----------------------
# First, define what "hitting the target" means. For example:
# Suppose current cumulative B sales are 1200, and the annual target is 1350.
# With 4 weeks left in the year, B needs to hit 37.5/week on average.
current_B_cumulative <- sum(df$B[df$date < as.Date("2024-01-01")])
target_total_B <- 1350
remaining_weeks <- 4
required_B_per_week <- (target_total_B - current_B_cumulative) / remaining_weeks

# Create the future path for B (matches target)
future_B <- rep(required_B_per_week, remaining_weeks)

# Generate future exogenous variables (Fourier terms + holidays)
future_dates <- seq.Date(as.Date("2024-01-01"), length.out = remaining_weeks, by = "week")
future_fourier <- fourier(ts_data, K = 2, h = remaining_weeks)
future_holidays <- data.frame(
  christmas = ifelse(week(future_dates) %in% c(51,52), 1, 0),
  thanksgiving = ifelse(week(future_dates) == 45, 1, 0)
)
future_exog <- cbind(future_fourier, future_holidays)
colnames(future_exog) <- colnames(exog_vars)

# Run conditional prediction: fix B's future path, predict A
conditional_pred <- predict(
  var_model,
  n.ahead = remaining_weeks,
  exogen = future_exog,
  condition = list(B = future_B)
)

# View the predicted A values
print(conditional_pred$fcst$A)

# Plot the forecast
plot(conditional_pred, main = "Conditional Forecast: Series A if Series B Hits Annual Target")

Critical Tips for Your Scenario #

• Small Sample Caution: Stick to lower lag orders (BIC will help with this) to avoid overfitting. If the model diagnostics look off, try reducing the lag count manually.
• Validate Scenarios: Test multiple versions of "B hitting target" (e.g., best-case vs. steady-case weekly B values) to see how sensitive A's forecast is.
• Alternative: VECM: If your series are cointegrated (use ca.jo() from the urca package to test), a VECM might perform better than VAR—it accounts for long-term equilibrium between A and B while modeling short-term changes.

内容的提问来源于stack exchange，提问作者dss