1) Problem 4.6

#import data
ex1 <- read_excel("GBweb_Row_Format.xlsx", sheet = "gRGDP", range = "A16:C153")
colnames(ex1) <- c('date', 'actual_growth', 'greenbook_growth_forecast')

Time Series Plots

#plot of realized values
actual_growth_ts <- ts(ex1$actual_growth, start = 1970, freq = 4)
plot(actual_growth_ts, main="Actual Growth", xlab = "Year", ylab = "Rate")

#plot of forecasts
greenbook_growth_forecast_ts <- ts(ex1$greenbook_growth_forecast, start = 1970, freq = 4)
plot(greenbook_growth_forecast_ts, main="Greenbook Growth Forecast", xlab = "Year", ylab = "Rate")

#plot of forecast errors
forecast_errors <- actual_growth_ts - greenbook_growth_forecast_ts
plot(forecast_errors, main="Forecast Errors", xlab = "Year", ylab = "Rate")

The realized GDP series and the Greenbook growth forecast exhibit similar cyclical patterns, repeating roughly every 15 years. In contrast, the time series of the forecast error displays stationarity of the first order, meaning it reverts to a constant mean over time. The forecast error spikes in 1985, 1995, and 2000 which correspond to instances where the Greenbook underestimated actual GDP growth.

Descriptive Statistics

summary(actual_growth_ts)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -10.900  -0.200   3.500   2.599   5.900   9.000
summary(greenbook_growth_forecast_ts)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -4.900   0.000   3.700   2.955   5.500   8.500
summary(forecast_errors)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -8.4000 -1.5000 -0.3000 -0.3555  0.7000  9.8000

The actual GDP growth series exhibits a wider range of values compared to the Greenbook forecasts, ranging from -10.9% to 9.0%. This suggests a higher degree of volatility in actual economic growth compared to the forecasts. The Greenbook growth forecasts show a slightly higher average (2.96%) compared to the actual growth (2.60%). The forecast errors have a mean close to zero (-0.36), indicating that the Greenbook forecasts are unbiased on average. However, the presence of negative minimum values suggests that the Greenbook has overestimated GDP growth in some periods.

ACF/PACF

par(mfrow=c(1,2))
#acf/pacf of realized values
acf(actual_growth_ts)
pacf(actual_growth_ts)

#acf/pacf of forecasts
acf(greenbook_growth_forecast_ts)
pacf(greenbook_growth_forecast_ts)

#acf/pacf of forecast errors
acf(forecast_errors)
pacf(forecast_errors)

par(mfrow=c(1,1))

The autocorrelation function (ACF) for GDP growth, Greenbook growth forecasts, and forecast errors all indicate dependence of present values on past values. However, this dependence weakens as the lag increases, with most autocorrelation falling within the significance level by the third lag.

2) Problem 4.8

#expected value of the forecast error
mean(forecast_errors)
## [1] -0.3554745

The expected value of the forecast errors is -0.3554745.

#run regression of forecast errors on several lags
reg <- dynlm(forecast_errors ~ forecast_errors + L(forecast_errors, 1) + L(forecast_errors, 2) + L(forecast_errors, 3))
## Warning in model.matrix.default(mt, mf, contrasts): the response appeared on the
## right-hand side and was dropped
## Warning in model.matrix.default(mt, mf, contrasts): problem with term 1 in
## model.matrix: no columns are assigned
summary(reg)
## 
## Time series regression with "ts" data:
## Start = 1970(4), End = 2004(1)
## 
## Call:
## dynlm(formula = forecast_errors ~ forecast_errors + L(forecast_errors, 
##     1) + L(forecast_errors, 2) + L(forecast_errors, 3))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.9385 -1.1059 -0.0343  0.8542  8.1279 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           -0.19979    0.20509  -0.974    0.332    
## L(forecast_errors, 1)  0.54642    0.08843   6.179  7.7e-09 ***
## L(forecast_errors, 2)  0.03297    0.10067   0.328    0.744    
## L(forecast_errors, 3) -0.06257    0.08946  -0.699    0.486    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.352 on 130 degrees of freedom
## Multiple R-squared:  0.2895, Adjusted R-squared:  0.2731 
## F-statistic: 17.65 on 3 and 130 DF,  p-value: 1.129e-09

The intercept is not statistically significant (p-value = 0.332 > 0.05), meaning the average forecast error is not statistically significant. The coefficient for L(forecast_errors, 1) (lagged error at period 1) is highly significant (p-value < 0.001) and positive (0.546). The coefficients for L(forecast_errors, 2) and L(forecast_errors, 3) (lagged errors at periods 2 and 3) are not statistically significant (p-value greater than 0.05)

#perform an f-test
f_test <- linearHypothesis(reg, c("L(forecast_errors, 1) = 0", "L(forecast_errors, 2) = 0", "L(forecast_errors, 3) = 0"), test = "F")
f_test

The F-statistic is 17.65457 and the p-value is 1.1289356e-09. Since the p-value is very low (essentially 0), we can reject the null hypothesis. Forecaset error has a strong correlation with its own past.

3) Problem 6.4

a) Theoretical autocorrelation for lag 1 \[ y_t = 0.7 - 2\epsilon_{t-1} + 1.35\epsilon_{t-2} + \epsilon_t \\ \rho_1 = \frac{\gamma_1}{\gamma_0} = \frac{\sigma^2(\theta_1 + \theta_1\theta_2^2)} {\sigma^2(1 + \theta_1^2 + \theta_2^2)} \\ = \frac{1(-2 + (-2)(1.35))} {1 (1+ (-2)^2)+(1.35)^2} \\ = \frac{4.77}{6.82} \\ = -0.689 \] b)

ma.sim <- arima.sim(model = list(ma=c(-2, 1.35)), n = 100)
acf(ma.sim)

pacf(ma.sim)

The autocorrelation seems to have an irregular cyclical pattern.

4) Problem 6.5-a

ma.sim<-arima.sim(model=list(ma=c(-2, 1.35)),n=100) 
quartz()
par(mfrow=(c(3,1)))
plot(ma.sim)
acf(ma.sim)
pacf(ma.sim)

The autocorrelation function exhibits values close to the theoretical correlation, but with a repeating pattern. Every four lags, the autocorrelation seems to return near zero. The only exception is the value for lag two, which is lower than the theoretical correlation at that lag.

4) Problem 6.5-b.

plot(forecast(ma.sim,h=1)) #One day

plot(forecast(ma.sim,h=2)) #Two day

plot(forecast(ma.sim,h=3)) #Three day

5) Problem 6.6

First AR Representation \[ \begin{split} y_t = 1.2 + 0.8\epsilon_{t-1} + \epsilon_t \\ y_t = 1.2 + (1 + 0.8L)\epsilon_t \\ 0.8 < 1 \\ \frac{1}{1+0.8L} = 1 - 0.8L + 0.8^2L^2-0.8^3L^3 + ... \\ \end{split} \]

Because 0.8 < 1, it is invertible. The autoregressive representation is:

\[ \begin{split} (y_t - 1.2) = 0.8(y_{t-1}-1.2) - 0.8^2(y_{t-2}-1.2) + 0.8^3(y_{t-3}-1.2)... + \epsilon_t\\ \end{split} \]

Second AR Representation

\[ \begin{split} y_t = 1.2 + 1.25\epsilon_{t-1} + \epsilon_t \\ |1.25| > 1 \\ \frac{1}{1+1.25L} = 0.8F * (1 - 0.8F + 0.8^2F^2 - 0.8^3F^3 + ...) \\ \end{split} \]

Because 1.25 > 1, it is non-invertible. The autoregressive representation is:

\[ \begin{split} (y_t - 1.2) * 0.8F(1 - 0.8F + 0.8^2F^2 - 0.8^3F^3 + ...) = \epsilon_t \end{split} \] The first representation is better because it includes how past errors contribute to the present value.

6) Problem 6.10

nvda <- read.csv("NVDA2.csv")
nvda_ts <- ts(nvda$AdjClose, start = 2019, end = 2024, frequency = 251)
plot(nvda_ts)

acf(nvda_ts, lag.max=50, main = "ACF of Nvidia Stock Prices")

Given the stochastic trend in the data, the MA process is likely a better fit compared to the AR process to capture the volatility, particularly the noticeable sharp rise in early 2020 followed by the subsequent sharp decline. I’ll compare an AR model with an MA model to confirm this.

ar1=arma(nvda_ts,order=c(1,0))
summary(ar1)
## 
## Call:
## arma(x = nvda_ts, order = c(1, 0))
## 
## Model:
## ARMA(1,0)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.743761 -0.043911 -0.006683  0.038476  0.598190 
## 
## Coefficient(s):
##            Estimate  Std. Error  t value Pr(>|t|)    
## ar1        0.995058    0.002637  377.409   <2e-16 ***
## intercept  0.010159    0.005750    1.767   0.0773 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Fit:
## sigma^2 estimated as 0.01264,  Conditional Sum-of-Squares = 15.85,  AIC = -1921.21
ar2=arma(nvda_ts,order=c(2,0))
summary(ar2)
## 
## Call:
## arma(x = nvda_ts, order = c(2, 0))
## 
## Model:
## ARMA(2,0)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.744271 -0.044028 -0.006686  0.038675  0.597914 
## 
## Coefficient(s):
##             Estimate  Std. Error  t value Pr(>|t|)    
## ar1        0.9959510   0.0282230   35.289   <2e-16 ***
## ar2       -0.0008814   0.0282069   -0.031   0.9751    
## intercept  0.0101096   0.0057574    1.756   0.0791 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Fit:
## sigma^2 estimated as 0.01265,  Conditional Sum-of-Squares = 15.85,  AIC = -1918.28
ar5=arma(nvda_ts,order=c(5,0))
summary(ar5)
## 
## Call:
## arma(x = nvda_ts, order = c(5, 0))
## 
## Model:
## ARMA(5,0)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.732180 -0.042576 -0.006815  0.038079  0.595544 
## 
## Coefficient(s):
##            Estimate  Std. Error  t value Pr(>|t|)    
## ar1        0.996110    0.028212   35.308   <2e-16 ***
## ar2       -0.037354    0.039827   -0.938   0.3483    
## ar3        0.035191    0.039828    0.884   0.3769    
## ar4        0.029733    0.039824    0.747   0.4553    
## ar5       -0.028616    0.028197   -1.015   0.3102    
## intercept  0.010190    0.005773    1.765   0.0776 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Fit:
## sigma^2 estimated as 0.01265,  Conditional Sum-of-Squares = 15.82,  AIC = -1912.2
quartz()

AR(5) model would be the best fit because it has the lowest AIC value.

plot(nvda_ts,xlab='Year', ylab="Nvdia", lwd=2)
grid()
lines(ar1$fitted.values,col="blue",lwd=2,lty=2)
lines(ar2$fitted.values,col="seagreen2",lwd=2,lty=2)
lines(ar5$fitted.values,col="red",lwd=2,lty=2)
legend("topright",legend=c("Data","AR(5)","AR(2)","AR(1)"),text.col=1:4,bty="n")

ar5=Arima(nvda_ts,order=c(5,0,0))
plot(forecast(ar5,h=20),xlim=c(2019,2024)) #20-days ahead

plot(forecast(ar5,h=100),xlim=c(2019,2024)) #100-days ahead

plot(forecast(ar5,h=200),xlim=c(2019,2024)) #200-days ahead

#MA(q) Model Fitting
ma1=arma(nvda_ts,order=c(0,1))
summary(ma1)
## 
## Call:
## arma(x = nvda_ts, order = c(0, 1))
## 
## Model:
## ARMA(0,1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.3966 -0.5145 -0.1285  0.4446  2.1682 
## 
## Coefficient(s):
##            Estimate  Std. Error  t value Pr(>|t|)    
## ma1        0.931970    0.007504   124.20   <2e-16 ***
## intercept  1.812244    0.035037    51.72   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Fit:
## sigma^2 estimated as 0.4156,  Conditional Sum-of-Squares = 521.19,  AIC = 2465.58
ma2=arma(nvda_ts,order=c(0,2)) 
summary(ma2)
## 
## Call:
## arma(x = nvda_ts, order = c(0, 2))
## 
## Model:
## ARMA(0,2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.40883 -0.30698 -0.08334  0.25745  1.72426 
## 
## Coefficient(s):
##            Estimate  Std. Error  t value Pr(>|t|)    
## ma1         1.47808     0.02145    68.92   <2e-16 ***
## ma2         0.80387     0.01458    55.15   <2e-16 ***
## intercept   1.79223     0.03747    47.83   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Fit:
## sigma^2 estimated as 0.1667,  Conditional Sum-of-Squares = 208.94,  AIC = 1320
ma5=arma(nvda_ts,order=c(0,5)) 
summary(ma5)
## 
## Call:
## arma(x = nvda_ts, order = c(0, 5))
## 
## Model:
## ARMA(0,5)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.34391 -0.13534 -0.03802  0.12243  1.05387 
## 
## Coefficient(s):
##            Estimate  Std. Error  t value Pr(>|t|)    
## ma1         1.67197     0.03165    52.82   <2e-16 ***
## ma2         1.97205     0.05029    39.21   <2e-16 ***
## ma3         1.73602     0.04465    38.88   <2e-16 ***
## ma4         1.14216     0.03304    34.57   <2e-16 ***
## ma5         0.45993     0.02468    18.64   <2e-16 ***
## intercept   1.70701     0.04708    36.26   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Fit:
## sigma^2 estimated as 0.04609,  Conditional Sum-of-Squares = 57.88,  AIC = -288.51

All the coefficients in the MA models are statistically significant. MA(5) model would be the best fit because it has the lowest AIC value. However, the lower AIC values in the MA models compared to the AR model indicate that the MA models better capture the data’s underlying structure, leaving less unexplained random shock.

#MA(q) fits
quartz()
plot(nvda_ts,xlab='Year', ylab="NVIDIA Stock Price", lwd=2)
grid()
lines(ma1$fitted.values,col="blue",lwd=1)
lines(ma2$fitted.values,col="seagreen2",lwd=1)
lines(ma5$fitted.values,col="red",lwd=1)
legend("topright",legend=c("Data","MA(5)","MA(2)","MA(1)"),text.col=1:4,bty="n")

#Forecast
ma5=Arima(nvda_ts,order=c(0,0,5))
plot(forecast(ma5,h=50)) #50 days ahead

plot(forecast(ma5,h=200)) #200 days ahead

7) Problem 6.2

a

library(fma)
autoplot(plastics) + ylab("Sales (Thousands)") + ggtitle("Monthly Sales of Product A for Plastics Manufacturer")

ggseasonplot(plastics)

ggsubseriesplot(plastics)

Sales of product A exhibit a seasonal trend, peaking in summer and showing an overall upward trend over five years.

b

#perform multiplicative decomposition
plastics_dcmp <- decompose(plastics, "multiplicative")

#plot multiplicative decomposition
autoplot(plastics_dcmp) + ggtitle("Sales of Product A for Plastics Manufacturer — Multiplicative Decomposition")

c

The multiplicative decomposition reinforces the graphical findings, highlighting the recurring summer highs in sales.

d

seasonal_plastics <- plastics_dcmp$seasonal
seasonally_adjusted_plastics <- plastics/seasonal_plastics

autoplot(plastics, series = "Data") + autolayer(seasonally_adjusted_plastics, series = "Seasonally Adjusted") + ggtitle("Sales of Product A for Plastics Manufacturer") + ylab("Sales (Thousands)")

#remove the trend and seasonality
trend_plastics <- plastics_dcmp$trend
detrend_plastics <- plastics/trend_plastics
seasonal_plastics <- plastics_dcmp$seasonal
seasonally_adjusted_plastics <- plastics/seasonal_plastics

#detrended and seasonally adjusted series
detrend_seas_adj_plastics <- plastics/plastics_dcmp$trend/plastics_dcmp$seasonal
autoplot(detrend_seas_adj_plastics, main = "Detrended and Seasonally Adjusted - Multiplicative")

#compute residuals
residuals_plastics <- plastics / (trend_plastics * seasonal_plastics)

# Plot ACF and PACF of residuals
tsdisplay(residuals_plastics, lag.max = 20, main = "ACF and PACF of Residuals - Multiplicative")

# Seasonally Adjust: Y / S
autoplot(seasonally_adjusted_plastics)

e

plastics2 <- plastics
plastics2[30] <- plastics2[30] + 500

plastics_dcmp2 <- decompose(plastics2, "multiplicative")

seasonal_plastics2 <- plastics_dcmp2$seasonal
seasonally_adjusted_plastics2 <- plastics2/seasonal_plastics2

autoplot(plastics2, series = "Data") + autolayer(seasonally_adjusted_plastics2, series = "Seasonally Adjusted") + ggtitle("Sales of Product A for Plastics Manufacturer (with outlier)") + ylab("Sales (Thousands)")

The outlier causes a significant upward spike in the data between 3 and 4.

plastics3 <- plastics
plastics3[55] <- plastics3[55] + 500

plastics_dcmp3 <- decompose(plastics3, "multiplicative")

seasonal_plastics3 <- plastics_dcmp3$seasonal
seasonally_adjusted_plastics3 <- plastics3/seasonal_plastics3

autoplot(plastics3, series = "Data") + autolayer(seasonally_adjusted_plastics3, series = "Seasonally Adjusted") + ggtitle("Sales of Product A for Plastics Manufacturer (with outlier)") + ylab("Sales (Thousands)")

The outlier’s placement at the end, rather than the middle, creates a spike in the data’s tail, not its center. This seems to have a greater impact on the overall upward trend than on the seasonal pattern, compared to a mid-series outlier.

8) Problem 6.6

a.

# STL decomposition with fixed seasonality
stl_brick_fixed_st <- stl(bricksq, 
                          s.window = "periodic",
                          robust = TRUE)

# STL decomposition with changing seasonality
stl_brick_changing_st <- stl(bricksq,
                             s.window = 5,
                             robust = TRUE)

# plot decomposed data
autoplot(stl_brick_fixed_st) +
  ggtitle("Brick production data decomposed by STL with fixed seasonality")

autoplot(stl_brick_changing_st) +
  ggtitle("Brick production data decomposed by STL with changing seasonality")

The plot indicates that the fixed seasonality exhibits an additive seasonal component, while the changing seasonality displays a multiplicative seasonal component. Additionally, the fixed seasonality shows a slightly noisier remainder compared to the changing seasonality.

b.

autoplot(bricksq, series = "Data") +
  autolayer(trendcycle(stl_brick_fixed_st),
            series = "Trend-cycle") +
  autolayer(seasadj(stl_brick_fixed_st),
            series = "Seasonally Adjusted Data") +
  ggtitle("Quarterly clay brick production in Australia",
          subtitle = "-decomposed by STL with fixed seasonality") +
  scale_color_manual(values = c("gray", "red", "blue"),
                     breaks = c("Data", "Trend-cycle", "Seasonally Adjusted Data"))

autoplot(bricksq, series = "Data") +
  autolayer(trendcycle(stl_brick_fixed_st),
            series = "Trend-cycle") +
  autolayer(seasadj(stl_brick_fixed_st),
            series = "Seasonally Adjusted Data") +
  ggtitle("Quarterly clay brick production in Australia",
          subtitle = "-decomposed by STL with changing seasonality") +
  scale_color_manual(values = c("gray", "red", "blue"),
                     breaks = c("Data", "Trend-cycle", "Seasonally Adjusted Data"))

c.

stl_brick_fixed_st %>% seasadj() %>% naive() %>% autoplot() + 
  ggtitle(label = "Naive forecast of seasonally adjusted brick data",
          subtitle = "after STL decomposition with fixed seasonality")

stl_brick_changing_st %>% seasadj() %>% naive() %>% autoplot() + 
  ggtitle(label = "Naive forecast of seasonally adjusted brick data",
          subtitle = "after STL decomposition with changing seasonality")

Comparing the two forecasts, we observe that the prediction intervals of seasonally adjusted data decomposed by STL with changing seasonality are slightly narrower than those with fixed seasonality. This is due to a reduction in the variance of the remainder component when the seasonality is allowed to change.

d. 

stlf_brick <- stlf(bricksq)
autoplot(stlf_brick)

### e.

checkresiduals(stlf_brick)

## 
##  Ljung-Box test
## 
## data:  Residuals from STL +  ETS(M,N,N)
## Q* = 44.704, df = 8, p-value = 4.186e-07
## 
## Model df: 0.   Total lags used: 8

The residuals appear to be correlated with each other, as indicated by the low p-value, suggesting a lack of independence.

f.

stlf_brick_robust <- stlf(bricksq, robust = TRUE)
autoplot(stlf_brick_robust)

checkresiduals(stlf_brick_robust)

## 
##  Ljung-Box test
## 
## data:  Residuals from STL +  ETS(M,N,N)
## Q* = 27.278, df = 8, p-value = 0.000633
## 
## Model df: 0.   Total lags used: 8

After applying a robust seasonal-trend decomposition, the autocorrelations have generally decreased. However, there are still some remaining high values.

g.

trainset_brick <- subset(bricksq, 
                        end = length(bricksq) - 8)
testset_brick <- subset(bricksq,
                        start = length(bricksq) - 7)

snaive_brick <- snaive(trainset_brick)
stlf_brick_part <- stlf(trainset_brick, robust = TRUE)
autoplot(bricksq, series = "Original data") +
  geom_line(size = 1) +
  autolayer(stlf_brick_part, PI = FALSE, size = 1,
            series = "stlf") +
  autolayer(snaive_brick, PI = FALSE, size = 1,
            series = "snaive") +
  scale_color_manual(values = c("gray50", "blue", "red"),
                     breaks = c("Original data", "stlf", "snaive")) +
  scale_x_continuous(limits = c(1990, 1994.5)) +
  scale_y_continuous(limits = c(300, 600)) +
  guides(colour = guide_legend(title = "Data")) +
  ggtitle("Forecast from stlf and snaive functions") +
  annotate(
    "rect",
    xmin=1992.75,xmax=1994.5,ymin=-Inf,ymax=Inf,
    fill="lightgreen",alpha = 0.3
    )
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Warning: Removed 136 rows containing missing values or values outside the scale range
## (`geom_line()`).

From the plot, it’s evident that the forecasts generated by the “stlf” function closely resemble the original data compared to those produced by the “snaive” function. Unlike the “snaive” function, “stlf” incorporates both trend and adaptable seasonality. Since the test set exhibits a trend with seasonality, “stlf” outperformed “snaive” in predicting near-future brick production amounts.

9) Problem 6.7

autoplot(writing)

Given the presence of both an upward trend and yearly seasonality in the data, the “rwdrift” method is more suitable than the “naive” method.

fit <- stlf(writing, s.window = "periodic", robust = TRUE, method = "rwdrift")

autoplot(fit, main = "Forecast without Box-Cox Transformation") + xlab("Month") + ylab("Sales")

fit2 <- stlf(writing, s.window = "periodic", robust = TRUE,lambda = BoxCox.lambda(writing))

autoplot(fit2, main = "Forecast with Box-Cox Transformation") + xlab("Month") + ylab("Sales")

Comparing the two forecasts, the Box-Cox transformation appears to yield superior performance with narrower prediction intervals. This enhancement in forecasting is attributed to the variability in seasonality observed over time in the original dataset.

10) Problem 6.8

str(fancy)
##  Time-Series [1:84] from 1987 to 1994: 1665 2398 2841 3547 3753 ...
head(fancy)
##          Jan     Feb     Mar     Apr     May     Jun
## 1987 1664.81 2397.53 2840.71 3547.29 3752.96 3714.74
autoplot(fancy)

Given the upward trend in the data, employing the rwdrift method for forecasting the non-seasonal component is preferable. Additionally, applying the Box-Cox transformation would be advantageous as it equalizes the variability in the data across the entire series.

stlf_fancy <- stlf(fancy,
                   s.window = 13,
                   robust = TRUE,
                   lambda = BoxCox.lambda(fancy),
                   method = "rwdrift")

autoplot(stlf_fancy)

The Box-Cox transformation significantly enlarges the prediction intervals. However, without this transformation, the forecasts become impractical.