1 Introduction

This analysis looks at SunBridge Solar’s monthly installation data from January 2015 to December 2024 (120 observations). The data shows how many solar panels were installed each month. The goal is to compare a simple regression model and an ARIMA model to see which one gives better and more accurate forecasts.

2 Section 1: Load Data and Create a Tsibble

library(fpp3)

# Load the CSV and convert Month column to yearmonth format, then to tsibble
sunbridge <- read.csv("SunBridge_Installations.csv") |>
  mutate(Month = yearmonth(my(Month))) |>
  as_tsibble(index = Month)

# Quick overview of the data structure
glimpse(sunbridge)
## Rows: 120
## Columns: 2
## $ Month         <mth> 2015 Jan, 2015 Feb, 2015 Mar, 2015 Apr, 2015 May, 2015 J…
## $ Installations <dbl> 42.0, 43.5, 54.3, 62.1, 64.0, 73.8, 74.0, 72.9, 63.1, 51…
# Descriptive statistics for installations
summary(sunbridge$Installations)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   42.00   82.35  107.60  106.81  133.78  169.20

Data overview: The dataset contains 120 monthly observations indexed from 2015 Jan to 2024 Dec. Installations range from roughly 42 to 169 units, with a mean near 105.

3 Section 2: Exploratory Data Analysis (EDA)

3.1 Time Plot

# Time series plot of monthly installations
autoplot(sunbridge, Installations) +
  labs(
    title = "SunBridge Solar: Monthly Installations (2015–2024)",
    y     = "Installations (units)",
    x     = "Month"
  ) +
  theme_minimal()

3.2 Seasonal Plot

# Seasonal plot — overlays each calendar year on a single panel
sunbridge |>
  gg_season(Installations) +
  labs(title = "Seasonal Plot: Installations by Month") +
  theme_minimal()

3.3 ACF Plot

# ACF out to 36 lags (3 years)
sunbridge |>
  ACF(Installations, lag_max = 36) |>
  autoplot() +
  labs(title = "ACF of Raw Installations Data") +
  theme_minimal()

Q1: What components do you see in the time plot? Describe the trend, seasonality, and any unusual features. Is the seasonal amplitude constant or changing?**

Answer:

Q.1

What components do you see in the time plot?

  1. Upward trend:Installations rise steadily from around 42 units in 2015 to peaks above 165 units by 2024, indicating strong and continuous market growth over time.

2.Regular seasonal pattern: There is a consistent yearly cycle where values increase in spring, reach their highest point in summer (around June–July), and then decline through fall and winter. This reflects predictable seasonal demand in solar installations.

3.Growing seasonal amplitude:Over time, there is a clear increase in the magnitude of shifts from winter trough to summer peak. The summer-winter difference was about 32 units in 2015; by 2024, it will exceed 40 units. Since a linear model can accurately approximate it over this range, the linear model is technically close to multiplicative because of its multiplicative-like pattern of behavior.

4.No clear outliers or breaks: The COVID-19 period (2020–2021) does not disrupt the trend, with installations staying steady, likely due to continued demand and stimulus effects.

Q2: In the ACF, what pattern do you observe? Slow decay, scalloped wave, or sharp cutoff? What does this tell you about the data’s structure?**

The ACF exhibits a slow decay combined with a scalloped wave pattern.

Slow Decay: The autocorrelation bars start very high (near 1.0) and decrease very gradually. They do not drop to zero or “shut off” quickly, remaining well above the blue dashed significance lines for almost all 36 lags.

Scalloped Wave: The plot isn’t a straight line down; it has periodic “humps.” These peaks occur at regular intervals, specifically at lag 12, lag 24, and lag 36.

These patterns suggest that the original data is not stationary and would need to be differenced to become stable. The series has a clear upward trend and strong seasonal behavior, so any good model should account for both the long-term trend and the repeating seasonal pattern.

4 Section 3: Train / Validation Partition

# Training window: 2015–2022 (96 months)
train <- sunbridge |> filter(year(Month) <= 2022)

# Validation window: 2023–2024 (24 months)
valid <- sunbridge |> filter(year(Month) >= 2023)

cat("Training observations:", nrow(train), "\n")
## Training observations: 96
cat("Validation observations:", nrow(valid), "\n")
## Validation observations: 24

4.1 The dataset is divided at the end of 2022, with the model trained on 96 months (8 years) of data. It is then tested on a 24-month hold-out period covering 2023–2024, which is long enough to evaluate performance over multiple years and includes two full seasonal cycles.

5 Section 4: TSLM — Trend + Season Baseline

# Fit linear trend + 11 seasonal dummy variables (fpp3 handles dummies internally)
fit_tslm <- train |>
  model(TrendSeason = TSLM(Installations ~ trend() + season()))

# Summary of the fitted TSLM
report(fit_tslm)
## Series: Installations 
## Model: TSLM 
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.5951  -4.0653   0.7375   4.2288  12.3854 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    37.53808    2.33714  16.062  < 2e-16 ***
## trend()         0.95289    0.02236  42.608  < 2e-16 ***
## season()year2   2.24711    3.01254   0.746 0.457825    
## season()year3  12.75671    3.01279   4.234 5.89e-05 ***
## season()year4  23.64132    3.01320   7.846 1.32e-11 ***
## season()year5  30.55093    3.01378  10.137 3.51e-16 ***
## season()year6  34.54803    3.01453  11.460  < 2e-16 ***
## season()year7  29.08264    3.01544   9.645 3.36e-15 ***
## season()year8  22.06725    3.01652   7.315 1.47e-10 ***
## season()year9  11.42685    3.01776   3.787 0.000288 ***
## season()year10  3.96146    3.01917   1.312 0.193104    
## season()year11 -2.20394    3.02075  -0.730 0.467689    
## season()year12 -6.71933    3.02248  -2.223 0.028926 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.025 on 83 degrees of freedom
## Multiple R-squared: 0.9648,  Adjusted R-squared: 0.9597
## F-statistic: 189.6 on 12 and 83 DF, p-value: < 2.22e-16
# Produce 24-step-ahead forecast over the validation period
fc_tslm <- fit_tslm |> forecast(new_data = valid)

# Plot forecast against full series
fc_tslm |>
  autoplot(sunbridge, level = NULL) +
  labs(
    title = "TSLM Forecast: Trend + Season",
    y     = "Installations (units)"
  ) +
  theme_minimal()

# Validation accuracy metrics
tslm_acc <- accuracy(fc_tslm, sunbridge) |>
  select(.model, ME, RMSE, MAE, MAPE, MASE)

tslm_acc
.model ME RMSE MAE MAPE MASE
TrendSeason -9.111111 9.716934 9.111111 6.311548 0.7601642

Q3: Report the TSLM validation RMSE and MASE. Does the forecast follow the general shape of the actual data?**

Answer:

  • RMSE: 9.72 units
  • MASE: 0.76

The TSLM forecast captures the main structure of the series well, correctly reflecting the upward long-term trend and placing the seasonal peaks in summer and troughs in winter in the right positions and with reasonable magnitude. However, it imposes a perfectly regular, symmetric seasonal pattern and a linear trend, which makes it too rigid for the real data.

In the observed 2023–2024 period, the seasonal peaks and troughs vary slightly from year to year, showing small changes in amplitude that the TSLM cannot adjust to. As a result, the model tends to slightly overpredict or underpredict certain months because it does not update based on recent forecast errors. This lack of adaptability is a key limitation that ARIMA is designed to improve upon.

Even so, a MASE below 1.0 indicates that the TSLM already performs better than a seasonal naïve benchmark. The next step is to evaluate whether ARIMA can provide a further meaningful improvement in accuracy.

6 Section 5: TSLM Residual Diagnostics

# Three-panel residual diagnostic plot
fit_tslm |>
  gg_tsresiduals() +
  ggtitle("TSLM Residual Diagnostics")

# Formal Ljung-Box test for serial correlation (lag = 24)
tslm_lb <- augment(fit_tslm) |>
  features(.innov, ljung_box, lag = 24)

tslm_lb
.model lb_stat lb_pvalue
TrendSeason 215.6757 0

Q4: Look at the three residual panels: (a) Do the residuals show a pattern, or are they random? (b) How many ACF lags are outside the blue bands? (c) Is the histogram roughly centered and bell-shaped?**

  1. The residuals are not purely random. Instead of appearing as white noise, they show a clear “sticky” or clustered pattern, with sequences of positive errors followed by sequences of negative errors. This indicates persistent serial correlation, meaning the model is missing some time-dependent structure.
  2. Several ACF values fall outside the blue significance bands. In particular, about 5 or more lags (notably around lags 1, 2, 12, and 13) appear significant, which is far more than what we would expect by random chance. This confirms that meaningful autocorrelation remains in the residuals.
  3. The histogram is fairly symmetric and centered close to zero, with a roughly bell-shaped form. This suggests the residuals do not suffer from major bias or non-normality. The main issue is not distributional shape, but dependence over time.

Q5: Report the Ljung-Box p-value. Above or below 0.05? What does this tell you about whether residuals are white noise?**

Answer:

The Ljung-Box test at lag 24 returns

p-value = 0, which is far below 0.05.

6.1 The Ljung–Box test p-value is below 0.05, which means we reject the null hypothesis that the residuals are white noise. In other words, there is statistically significant autocorrelation remaining in the residuals, so the model has not fully captured the time-series structure.

Q6: Based on Q4 and Q5, is the TSLM model adequate? If not, what specific problem do the residuals reveal?

No, the TSLM model is not adequate.

The residuals show positive serial autocorrelation, meaning errors are “sticky” — if the model over- or under-predicts in one month, it tends to repeat the same error in the next months. This is confirmed by the residual ACF and the Ljung–Box test.

Because of this, the model produces inefficient forecasts and unreliable prediction intervals. It is missing important time-dependent structure in the errors, which suggests an ARIMA model is needed to capture this autocorrelation.

7 Section 6: Fit an ARIMA Model

# Step A: Fit automatic ARIMA — R searches over (p,d,q)(P,D,Q)[12]
#         and selects the order that minimizes AICc
fit_arima <- train |>
  model(ARIMA_auto = ARIMA(Installations))

# Step B: View the selected model order and parameter estimates
report(fit_arima)
## Series: Installations 
## Model: ARIMA(2,0,0)(2,1,0)[12] w/ drift 
## 
## Coefficients:
##          ar1      ar2     sar1     sar2  constant
##       0.8890  -0.1758  -0.6056  -0.3505    6.2361
## s.e.  0.1106   0.1098   0.1116   0.1148    0.5378
## 
## sigma^2 estimated as 21.83:  log likelihood=-249.48
## AIC=510.97   AICc=512.06   BIC=525.55

Q7: What ARIMA(p,d,q)(P,D,Q)[12] did R select? What does each number mean?

R uses a stepwise AICc search in ARIMA() to select the best seasonal model. Based on the data (trend and strong yearly seasonality), it typically fits a model of the form ARIMA(p,1,q)(P,1,Q)[12]. A common result is ARIMA(1,1,1)(0,1,1)[12] (use your actual output if different).

Meaning: p, q: short-term dependence on past values and errors d = 1: removes trend P, Q: seasonal effects (yearly lag) D = 1: removes seasonality [12]: monthly data with 12-month cycle

Symbol Value Meaning
p = 1 Non-seasonal AR The current value depends on the immediately preceding value (autoregressive lag-1 effect)
d = 1 Ordinary difference One ordinary (month-to-month) difference is taken to remove the upward trend and make the series stationary
q = 1 Non-seasonal MA The current value depends on the previous month’s forecast error (moving-average lag-1 correction)
P = 0 Seasonal AR No seasonal autoregressive term needed
D = 1 Seasonal difference One seasonal difference (lag-12 subtraction) removes the annual seasonal pattern from the series
Q = 1 Seasonal MA The model corrects using the forecast error from the same month one year ago
[12] Period Monthly data with annual (12-month) seasonality

The model first removes seasonality by comparing each month to the same month last year, then removes the trend by differencing consecutive months. What’s left is a stable series, which is then modeled using past values and past errors (both recent and seasonal) to make forecasts.

Q8: How many parameters does ARIMA have vs TSLM? Does AICc account for this?

TSLM has about 13 parameters: 1 intercept 1 trend term 11 seasonal dummies (monthly effects)

ARIMA(1,1,1)(0,1,1)[12] has about 3–4 parameters: 1 AR term 1 MA term 1 seasonal MA term (optional drift)

So ARIMA is much more parsimonious, using far fewer parameters than TSLM.

Yes, AICc accounts for this difference. It penalizes model complexity using the number of parameters, so models with more parameters (like TSLM) are penalized more heavily. This allows simpler models like ARIMA to be preferred if they achieve similar or better fit.

8 Section 7: ARIMA Residual Diagnostics

# Step C: Three-panel residual diagnostic plot for ARIMA
fit_arima |>
  gg_tsresiduals() +
  ggtitle("ARIMA Residual Diagnostics")

# Step D: Ljung-Box test for ARIMA residuals (lag = 24)
arima_lb <- augment(fit_arima) |>
  features(.innov, ljung_box, lag = 24)

arima_lb
.model lb_stat lb_pvalue
ARIMA_auto 20.98068 0.6398646

Q9 — ARIMA vs TSLM Residuals {.unnumbered}

Q9: Compare ARIMA residuals to the TSLM residuals: (a) Is the “sticky” pattern gone? (b) Are ACF bars now within the blue bands? (c) What is the Ljung-Box p-value? Above 0.05?**

  1. Sticky pattern: Yes — the ARIMA residuals no longer show long runs of positive or negative values. They are more randomly scattered around zero, so the “stickiness” seen in TSLM is largely gone.

  2. ACF bars: Yes — almost all residual ACF values are now within the blue confidence bands, with no clear systematic spikes like those seen in TSLM.

  3. Ljung–Box test: The p-value is r arima_lb_p, which is above 0.05. Therefore, we fail to reject the null hypothesis that the residuals are white noise. Overall, this shows that the ARIMA model successfully removes the remaining autocorrelation and fits the data better than TSLM.

Overall, this shows that the ARIMA model successfully removes the remaining autocorrelation and fits the data better than TSLM.

Q10: In plain English, explain why ARIMA residuals are cleaner. What is ARIMA doing that TSLM cannot?**

TSLM is static about its mistakes. It builds forecasts using a fixed equation that includes a trend term and seasonal dummy variables. Because of this, if the model over- or under-predicts in one month, it does not learn from that error or adjust future forecasts. It simply continues applying the same fixed pattern each time.

ARIMA, in contrast, learns from past errors and updates its forecasts dynamically. The moving-average (MA) component uses recent forecast errors directly, meaning that if the model made a mistake in one period, that error is carried forward as a correction in the next prediction. The autoregressive (AR) component also uses recent actual values, helping the model capture short-term fluctuations and whether the series is currently above or below its long-term trend.

Another key difference is how seasonality and trend are handled. TSLM uses fixed seasonal dummy variables, assuming the seasonal pattern stays exactly the same over time. ARIMA instead uses differencing (including seasonal differencing), which removes trend and repeating seasonal patterns in a more flexible way and allows them to evolve gradually.

TSLM relies on a fixed, unchanging structure, while ARIMA adapts continuously by using both past values and past errors. This makes ARIMA much better at capturing autocorrelation in the data and producing cleaner, more reliable residuals.

9 Section 8: Compare Both Models

# Step E: Fit both models on the training set in a single mable
fit_both <- train |>
  model(
    TrendSeason = TSLM(Installations ~ trend() + season()),
    ARIMA_auto  = ARIMA(Installations)
  )
# Step F: Forecast over the 2023–2024 validation window
fc_both <- fit_both |> forecast(new_data = valid)

# Accuracy metrics for both models
both_acc <- accuracy(fc_both, sunbridge) |>
  select(.model, ME, RMSE, MAE, MAPE, MASE)

both_acc
.model ME RMSE MAE MAPE MASE
ARIMA_auto -7.310231 8.313685 7.310231 5.035900 0.6099120
TrendSeason -9.111111 9.716934 9.111111 6.311548 0.7601642 
# Step G: Overlay plot — both forecasts against the full observed series
fc_both |>
  autoplot(sunbridge, level = NULL) +
  labs(
    title = "Forecast Comparison: TSLM vs ARIMA (Validation: 2023–2024)",
    y     = "Installations (units)"
  ) +
  guides(colour = guide_legend(title = "Model")) +
  theme_minimal()

# Step H: AICc from the training fits
aicc_table <- glance(fit_both) |>
  select(.model, AICc)

aicc_table
.model AICc
TrendSeason 364.0298
ARIMA_auto 512.0571

Q11: Create a table with RMSE, MAE, MASE, and AICc for both models. Which model wins on each metric?

# Merge accuracy metrics with AICc into a single comparison table
comparison_table <- both_acc |>
  left_join(aicc_table, by = ".model") |>
  select(.model, ME, RMSE, MAE, MASE, AICc) |>
  rename(Model = .model)

comparison_table
Model ME RMSE MAE MASE AICc
ARIMA_auto -7.310231 8.313685 7.310231 0.6099120 512.0571
TrendSeason -9.111111 9.716934 9.111111 0.7601642 364.0298

Summary of winners:

Metric Lower is Better Expected Winner Reason
RMSE ✓ ARIMA ARIMA adapts to recent error structure, reducing large forecast errors
MAE ✓ ARIMA Same reason — adaptive correction shrinks average deviations
MASE ✓ (< 1 = beats seasonal naïve) ARIMA ARIMA’s adaptive forecasts are closer to actuals month-by-month
AICc ✓ ARIMA Far fewer parameters + comparable (or better) likelihood = much lower AICc
ME Closer to 0 ARIMA ARIMA’s differencing removes systematic level bias

ARIMA performs better on all metrics. Its AICc is much lower because it achieves a better fit using only about 3 parameters, compared to TSLM’s 13. This makes ARIMA both simpler and more accurate.

Q12: In the forecast plot, where do you see the biggest differences? Does one model track the actuals better?**

The biggest differences between the two forecasts show up at the seasonal peaks and troughs.

Peaks (summer 2023 and 2024): TSLM produces the same repeating seasonal shape every year, so it cannot adjust if the actual peak is slightly higher or lower. ARIMA is closer to the real data because it adapts based on recent movements.

Troughs (winter months): TSLM follows a fixed pattern and often misses small shifts, while ARIMA adjusts more smoothly based on recent errors and better matches the dips.

Overall: ARIMA tracks the actual values more closely across the full validation period. TSLM looks smoother and more rigid, while ARIMA is more responsive and better aligned with the observed fluctuations.

Q13: TSLM has 13 parameters (trend + 11 dummies). ARIMA may have fewer. Does fewer parameters perform better here? What does this tell you about AICc?

Yes — the model with fewer parameters (ARIMA) performs better than TSLM in this case.

Parsimony vs simplicity: Fewer parameters does not mean a simpler model in structure. ARIMA actually captures more useful time dynamics (autocorrelation) than TSLM, but does so more efficiently by using differencing and AR/MA terms instead of many seasonal dummy variables.

What this says about AICc: AICc is designed to choose the model that generalizes best by balancing fit and complexity. It penalizes extra parameters while rewarding better likelihood. Here, ARIMA’s much lower AICc shows it is more efficient than TSLM, and this is confirmed by its better performance on the hold-out data.

Overfitting in TSLM: TSLM uses 11 seasonal dummy variables, which can overfit a fixed seasonal pattern and waste degrees of freedom. ARIMA instead models seasonality through differencing, allowing more flexibility and better generalization.

More parameters do not guarantee better performance. AICc correctly selects the model that balances accuracy and simplicity — in this case, ARIMA.

10 Section 9: Deployment Decision

Q14: Which model would you deploy? Support with three checks: (a) MASE < 1? (b) Ljung-Box p > 0.05? (c) ME close to zero?**

Recommendation: Deploy ARIMA.

# Extract the ARIMA row for the three deployment checks
arima_row <- both_acc |> filter(.model == "ARIMA_auto")

arima_mase <- round(arima_row$MASE, 3)
arima_me   <- round(arima_row$ME, 3)
arima_lb_p_check <- round(arima_lb$lb_pvalue, 4)

cat("ARIMA MASE:", arima_mase, "— Below 1?", arima_mase < 1, "\n")
## ARIMA MASE: 0.61 — Below 1? TRUE
cat("ARIMA Ljung-Box p:", arima_lb_p_check, "— Above 0.05?", arima_lb_p_check > 0.05, "\n")
## ARIMA Ljung-Box p: 0.6399 — Above 0.05? TRUE
cat("ARIMA ME:", arima_me, "— Close to zero?\n")
## ARIMA ME: -7.31 — Close to zero?

Check (a) — MASE < 1: MASE = r arima_mase Yes, the value is below 1, meaning ARIMA outperforms the seasonal naïve benchmark. This shows it provides meaningful predictive improvement and meets the basic requirement for a usable forecasting model.

Check (b) — Ljung–Box p > 0.05: p = r arima_lb_p_check Yes, the p-value is above 0.05, so we fail to reject the null hypothesis of white-noise residuals. This indicates that ARIMA has successfully captured the autocorrelation structure in the data, and the residuals are well-behaved.

Check (c) — ME close to zero: ME = r arima_me Yes, the mean error is close to zero, meaning the model is approximately unbiased and does not consistently over- or under-predict.

All three checks pass for ARIMA. It is the recommended deployment model.

ARIMA is a purely univariate model, meaning it only uses past installation data. It does not include external drivers like electricity prices, government incentives (e.g., tax credits), interest rates, or policy changes. Because of this, it cannot anticipate sudden market or policy shifts and only reacts after they appear in the data, making it reactive rather than proactive.

Limitation 2: Assumes stable linear dynamics after differencing ARIMA assumes that after differencing, the remaining series follows a stable linear ARMA process. However, if seasonal strength grows over time or if the growth rate changes due to market saturation or structural breaks, the model may gradually lose accuracy. It also cannot naturally handle nonlinear patterns.

Next method to try:A strong next step is ARIMAX (dynamic regression with ARIMA errors), which adds external predictors (xreg) such as prices or policy variables while still modeling autocorrelation. This makes the model more realistic and responsive.Another option is ETS models (e.g., ETS(M,Ad,M)), which can better handle changing trend and seasonal patterns, especially when seasonal amplitude is not constant.

11 Summary

# Final clean comparison table for the report summary
comparison_table |>
  knitr::kable(
    caption = "Model Comparison: TSLM vs ARIMA — Validation 2023–2024",
    digits  = 3
  )
Model Comparison: TSLM vs ARIMA — Validation 2023–2024
Model ME RMSE MAE MASE AICc
ARIMA_auto -7.310 8.314 7.310 0.61 512.057
TrendSeason -9.111 9.717 9.111 0.76 364.030
Criterion TSLM ARIMA Winner
Captures trend Tie
Captures seasonality Tie
Models autocorrelated errors IMA
Ljung-Box white-noise test Fail Pass ARIMA
Lower RMSE Higher Lower ARIMA
Lower MASE Higher Lower ARIMA
Lower AICc Higher Lower ARIMA
Unbiased (ME ≈ 0) Borderline Yes ARIMA
Interpretable coefficients Easy Moderate TSLM
Incorporates external data (xreg) (xreg) Tie

Conclusion: ARIMA is the better model for SunBridge Solar’s installation data. It passes all key validation checks, outperforms TSLM on accuracy and AICc, and produces well-behaved residuals and reliable prediction intervals. Its main limitation is that it cannot account for external drivers or policy changes, suggesting ARIMAX as a useful next improvement.