Forecasting, Inverse Transformation, and Bias Correction

Chapter 5: Lesson 5

Learning Outcomes

Apply logarithmic transformations to time series

• Apply a log-transformation to a multiplicative time series

Apply the bias correction factor for inverse transformations

• State the bias correction procedure for log-transform estimates
• Explain when to use the bias correction factor
• Use the bias correction factor for a log-transform model
• Forecast using the inverse-transform and bias correction of a log-transformed model

Preparation

• Read Sections 5.9-5.11

Learning Journal Exchange (10 min)

• Review another student’s journal

• What would you add to your learning journal after reading another student’s?

• What would you recommend the other student add to their learning journal?

• Sign the Learning Journal review sheet for your peer

Class Activity: Anti-Log Transformation and Bias Correction on Simulated Data (10 min)

Forecasts for a Simulated Time Series

In the previous lesson, we simulated a time series with an exponential trend. In this lesson, we will forecast future values based on our model.

Figure 1 shows the simulated time series and the time series after the natural logarithm is applied.

set.seed(12345)

n_years <- 9 # Number of years to simulate
n_months <- n_years * 12 # Number of months
sigma <- .05 # Standard deviation of random term

z_t <- rnorm(n = n_months, mean = 0, sd = sigma)

dates_seq <- seq(floor_date(now(), unit = "year"), length.out=n_months + 1, by="-1 month") |>
  floor_date(unit = "month") |> sort() |> head(n_months)

sim_ts <- tibble(
  t = 1:n_months,
  dates = dates_seq,
  random = arima.sim(model=list(ar=c(.5,0.2)), n = n_months, sd = 0.02),
  x_t = exp(2 + 0.015 * t +
              0.03 * sin(2 * pi * 1 * t / 12) + 0.04 * cos(2 * pi * 1 * t / 12) + 
              0.05 * sin(2 * pi * 2 * t / 12) + 0.03 * cos(2 * pi * 2 * t / 12) +
              0.01 * sin(2 * pi * 3 * t / 12) + 0.005 * cos(2 * pi * 3 * t / 12) +
              random
  )
) |>
  mutate(
    cos1 = cos(2 * pi * 1 * t / 12),
    cos2 = cos(2 * pi * 2 * t / 12),
    cos3 = cos(2 * pi * 3 * t / 12),
    cos4 = cos(2 * pi * 4 * t / 12),
    cos5 = cos(2 * pi * 5 * t / 12),
    cos6 = cos(2 * pi * 6 * t / 12),
    sin1 = sin(2 * pi * 1 * t / 12),
    sin2 = sin(2 * pi * 2 * t / 12),
    sin3 = sin(2 * pi * 3 * t / 12),
    sin4 = sin(2 * pi * 4 * t / 12),
    sin5 = sin(2 * pi * 5 * t / 12),
    sin6 = sin(2 * pi * 6 * t / 12)) |>
  mutate(std_t = (t - mean(t)) / sd(t)) |>
  as_tsibble(index = dates)

sim_plot_raw <- sim_ts |>
  autoplot(.vars = x_t) +
  labs(
    x = "Month",
    y = "x_t",
    title = "Simulated Time Series"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5)
  )

sim_plot_log <- sim_ts |>
  autoplot(.vars = log(x_t)) +
  labs(
    x = "Month",
    y = "log(x_t)",
    title = "Logarithm of Simulated Time Series"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5)
  )

sim_plot_raw | sim_plot_log

Figure 1: Time plot of the time series (left) and the natural logarithm of the time series (right)

We can use the forecast() function to predict future values of this time series. Table 1 displays the output of the forecast() command. Note that the column labeled x_t (i.e. \(x_t\)), representing the time series is populated with information tied to a normal distribution. The mean and standard deviation specified are the estimated parameters for the distribution of the predicted values of \(\log(x_t)\). If you raise \(e\) to the power of the mean, you get the values in the .mean column.

# Fit model (OLS)
sim_reduced_linear_lm1 <- sim_ts |>
  model(sim_reduced_linear1 = TSLM(log(x_t) ~ std_t + 
        sin1 + cos1 + sin2 + cos2 + sin3 + cos3))

# Compute forecast
n_years_forecast <- 5
n_months_forecast <- 12 * n_years_forecast

new_dat <- tibble(t = n_months:(n_months + n_months_forecast )) |>
  mutate(
    dates = seq(max(dates_seq), length.out=n_months_forecast + 1, by="1 month") 
  ) |>
  mutate(
    std_t = (t - mean(pull(sim_ts, t))) / sd(pull(sim_ts, t)),
    sin1 = sin(2 * pi * 1 * t / 12),
    cos1 = cos(2 * pi * 1 * t / 12),
    sin2 = sin(2 * pi * 2 * t / 12),
    cos2 = cos(2 * pi * 2 * t / 12),
    sin3 = sin(2 * pi * 3 * t / 12),
    cos3 = cos(2 * pi * 3 * t / 12),
    sin4 = sin(2 * pi * 4 * t / 12),
    cos4 = cos(2 * pi * 4 * t / 12),
    sin5 = sin(2 * pi * 5 * t / 12),
    cos5 = cos(2 * pi * 5 * t / 12),
    cos6 = cos(2 * pi * 6 * t / 12)
  ) |>
  as_tsibble(index = dates)

sim_reduced_linear_lm1 |> 
  forecast(new_data = new_dat)

Table 1: Output of the forecast() command for the simulated time series

.model	dates	x_t	.mean	t	std_t	sin1	cos1	sin2	cos2	sin3	cos3
sim_reduced_linear1	2023-12-01	t(N(3.7, 0.00055))	40.714	108	1.708	0	1	0	1	0	1	...
sim_reduced_linear1	2024-01-01	t(N(3.8, 0.00056))	43.133	109	1.74	0.5	0.866	0.866	0.5	1	0	...
sim_reduced_linear1	2024-02-01	t(N(3.7, 0.00056))	41.625	110	1.772	0.866	0.5	0.866	-0.5	0	-1	...
sim_reduced_linear1	2024-03-01	t(N(3.7, 0.00056))	39.159	111	1.804	1	0	0	-1	-1	0	...
⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮
sim_reduced_linear1	2028-11-01	t(N(4.5, 6e-04))	90.307	167	3.592	-0.5	0.866	-0.866	0.5	-1	0	...
sim_reduced_linear1	2028-12-01	t(N(4.6, 6e-04))	101.132	168	3.624	0	1	0	1	0	1	...

Figure 2 illustrates the forecasted values for the time series.

sim_forecast_plot_regular <- sim_reduced_linear_lm1 |> 
  forecast(new_data = new_dat) |>
  autoplot(sim_ts, level = 95) +
  labs(
    x = "Month",
    y = "x_t",
    title = "Simulated Time Series"
  ) +
  theme_minimal() +
  theme(legend.position = "inset") +
  theme(
    plot.title = element_text(hjust = 0.5)
  )

sim_forecast_plot_logged <- sim_reduced_linear_lm1 |> 
  forecast(new_data = new_dat) |>
  autoplot(sim_ts, level = 95) +
  scale_y_continuous(trans = "log", labels = trans_format("log")) +
  labs(
    x = "Month",
    y = "log(x_t)",
    title = "Logarithm of Simulated Time Series"
  ) +
  theme_minimal() +
  theme(legend.position = "inset")  +
  theme(
    plot.title = element_text(hjust = 0.5)
  )

sim_forecast_plot_regular

Figure 2: Forecasted values of the time series with 95% confidence bands

Bias Correction

The forecasts presented above were computed by raising \(e\) to the power of the predicted log-values. Unfortunately, this introduces bias in the forecasted means. This bias tends to be large if the regression model does not fit the data closely.

The textbook points out that the bias correction should only be applied for means, not for simulated values. This means that if you are simulating transformed values, and you apply the inverse of your original transformation, the resulting values are appropriate.

When we apply the inverse transform to the residual series, we introduce a bias. We can account for this bias applying one of two adjustments to our mean values. The theory behind this transformations is alluded to in the textbook, but is not essential.

There are two common patterns observed in the residual series: (1) Gaussian white noise or (2) Negatively-skewed values. Note that negatively-skewed values are the same as left-skewed values.

We can use the skewness statistic to assess the shape of the residual series. When the skewness is less than -1 or greater than 1, we say that the distribution is highly skewed. For skewness values between -1 and -0.5 or between 0.5 and 1, we say there is moderate skewness. If skewness lies between -0.5 and 0.5, the distribution is considered roughly symmetric.

Log-Normal Correction

Normally-Distributed Residual Series

If the residual series follows a normal distribution, we multiply the means of the forecasted values \(\hat x_t\) by the factor \(e^{\frac{1}{2} \sigma^2}\):

\[ \hat x_t' = e^{\frac{1}{2} \sigma^2} \cdot \hat x_t \]

where \(\left\{ \hat x_t: t = 1, \ldots, n \right\}\) gives the values of the forecasted series, and \(\left\{ \hat x_t': t = 1, \ldots, n \right\}\) is the adjusted forecasted values.

Emperical Correction

Negatively-Skewed Residual Series

If the residual series demonstrates negative skewness, then we can adjust the forecasts \(\left\{ \hat x_t \right\}\) as follows:

\[ \hat x_t' = e^{\widehat{\log x_t}} \sum_{t=1}^{n} \frac{e^{z_t}}{n} \]

where \(\left\{ \widehat{\log x_t}: t = 1, \ldots, n \right\}\) is the forecasted series obtained by fitting the log-regression model.

\(\left\{ z_t \right\}\) is the residual series from this fitted model in the log-transformed units.

The code given below can be used to compute the corrected mean values for the simulated data.

sim_model_values <- sim_reduced_linear_lm1 |>
  glance()

sim_model_check <- sim_model_values |>
  mutate(
    sigma = sqrt(sigma2),
    lognorm_cf = exp((1/2) * sigma2),
    empirical_cf = sim_reduced_linear_lm1 |>
      residuals() |>
      pull(.resid) |>
      exp() |>
      mean()) |>
  select(.model, r_squared, sigma2, sigma, lognorm_cf, empirical_cf)

sim_pred <- sim_reduced_linear_lm1 |> 
  forecast(new_data = new_dat) |>
  mutate(.mean_correction = .mean * sim_model_check$empirical_cf) |>
  select(t, x_t, .mean, .mean_correction)

sim_model_check

# A tibble: 1 × 6
  .model              r_squared   sigma2  sigma lognorm_cf empirical_cf
  <chr>                   <dbl>    <dbl>  <dbl>      <dbl>        <dbl>
1 sim_reduced_linear1     0.998 0.000507 0.0225       1.00         1.00

From this, we observe that for the simulated data, \(R^2 = 0.998\). This indicates that the model explains a high proportion of the variation in the data. The log-normal adjustment is \(1.00025\), and the emperical adjustment is \(1.00023\). Both of these values are extremely close to 1, so they will have a negligible impact on the predicted values.

This result does not generalize. In other situations, there can be a substantial effect of this bias on the predicted means.

Histogram of residuals

Figure 3 gives a histogram of the residuals and compute the skewness of the residual series.

sim_resid_df <- sim_reduced_linear_lm1 |> 
  residuals() |> 
  as_tibble() |> 
  dplyr::select(.resid) |>
  rename(x = .resid) 
  
sim_resid_df |>
  mutate(density = dnorm(x, mean(sim_resid_df$x), sd(sim_resid_df$x))) |>
  ggplot(aes(x = x)) +
    geom_histogram(aes(y = after_stat(density)),
        color = "white", fill = "#56B4E9", binwidth = 0.01) +
    geom_line(aes(x = x, y = density)) +
    theme_bw() +
    labs(
      x = "Values",
      y = "Frequency",
      title = "Histogram of Residuals from the Reduced Linear Model"
    ) +
    theme(
      plot.title = element_text(hjust = 0.5)
    )

Figure 3: Histogram of the values in the residual series based on the model with a linear trend and seasonal Fourier terms where i≤3

We can use the command skewness(sim_resid_df$x) to compute the skewness of these residuals: -0.135. This number is close to zero (specifically between -0.5 and 0.5,) so we conclude that the residual series is approximately normally distributed. We can apply the log-normal correction to our mean forecast values.

Class Activity: Apple Revenue (10 min)

We take another look at the quarterly revenue reported by Apple Inc. from Q1 of 2005 through Q1 of 2012

Visualizing the Time Series

Figure 4 gives the time plot illustrating the quarterly revenue reported by Apple from the first quarter of 2005 through the first quarter of 2012.

apple_raw <- rio::import("https://byuistats.github.io/timeseries/data/apple_revenue.csv") |>
  mutate(dates = round_date(mdy(date), unit = "quarter")) |>
  arrange(dates)

apple_ts <- apple_raw |>
  filter(dates <= my("Jan 2012")) |>
  dplyr::select(dates, revenue_billions) |>
  mutate(t = 1:n()) |>
  mutate(std_t = (t - mean(t)) / sd(t)) |>
  mutate(
    sin1 = sin(2 * pi * 1 * t / 4),
    cos1 = cos(2 * pi * 1 * t / 4),
    cos2 = cos(2 * pi * 2 * t / 4)
  ) |>
  as_tsibble(index = dates)

apple_plot_regular <- apple_ts |>
  autoplot(.vars = revenue_billions) +
  stat_smooth(method = "lm", 
              formula = y ~ x, 
              geom = "smooth",
              se = FALSE,
              color = "#E69F00",
              linetype = "dotted") +
    labs(
      x = "Quarter",
      y = "Revenue (Billions USD)",
      title = "Apple Revenue (in Billions USD)"
    ) +
    theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5))

apple_plot_transformed <- apple_ts |>
  autoplot(.vars = log(revenue_billions)) +
  stat_smooth(method = "lm", 
              formula = y ~ x, 
              geom = "smooth",
              se = FALSE,
              color = "#E69F00",
              linetype = "dotted") +
    labs(
      x = "Quarter",
      y = "Logarithm of Revenue",
      title = "Logarithm of Apple Revenue"
    ) +
    theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5))

apple_plot_regular | apple_plot_transformed

Figure 4: Apple quarterly revenue figures (in billions of U.S. dollars) from Q1 of 2005 to Q1 of 2012; the figure on the left presents the revenue in dollars and the figure on the right gives the logarithm of the quarterly revenue; a simple linear regression line is given for reference

Finding a Suitable Model

We start by fitting a cubic trend to the logarithm of the quarterly revenues. The full model is fitted here:

# Cubic model with standardized time variable

apple_full_cubic_lm <- apple_ts |>
  model(apple_full_cubic = TSLM(log(revenue_billions) ~ std_t + I(std_t^2) + I(std_t^3) +
    sin1 + cos1 + cos2 )) # Note sin2 is omitted
apple_full_cubic_lm |>
  tidy() |>
  mutate(sig = p.value < 0.05)

# A tibble: 7 × 7
  .model           term        estimate std.error statistic  p.value sig  
  <chr>            <chr>          <dbl>     <dbl>     <dbl>    <dbl> <lgl>
1 apple_full_cubic (Intercept)   1.83      0.0572    32.1   5.67e-20 TRUE 
2 apple_full_cubic std_t         1.01      0.0972    10.4   6.41e-10 TRUE 
3 apple_full_cubic I(std_t^2)   -0.0158    0.0445    -0.355 7.26e- 1 FALSE
4 apple_full_cubic I(std_t^3)    0.0866    0.0516     1.68  1.07e- 1 FALSE
5 apple_full_cubic sin1          0.217     0.0533     4.08  4.96e- 4 TRUE 
6 apple_full_cubic cos1          0.0839    0.0552     1.52  1.43e- 1 FALSE
7 apple_full_cubic cos2         -0.0981    0.0382    -2.57  1.74e- 2 TRUE 

The quadratic and cubic trend terms are not statistically signficant. We now eliminate the cubic term and fit a full model with a quadratic trend.

apple_full_quad_lm <- apple_ts |>
  model(apple_full_quad = TSLM(log(revenue_billions) ~ std_t + I(std_t^2) + 
    sin1 + cos1 + cos2 )) # Note sin2 is omitted
apple_full_quad_lm |>
  tidy() |>
  mutate(sig = p.value < 0.05) 

# A tibble: 6 × 7
  .model          term        estimate std.error statistic  p.value sig  
  <chr>           <chr>          <dbl>     <dbl>     <dbl>    <dbl> <lgl>
1 apple_full_quad (Intercept)   1.83      0.0594    30.9   3.12e-20 TRUE 
2 apple_full_quad std_t         1.16      0.0403    28.7   1.64e-19 TRUE 
3 apple_full_quad I(std_t^2)   -0.0158    0.0462    -0.341 7.36e- 1 FALSE
4 apple_full_quad sin1          0.217     0.0554     3.93  6.73e- 4 TRUE 
5 apple_full_quad cos1          0.0934    0.0570     1.64  1.15e- 1 FALSE
6 apple_full_quad cos2         -0.0981    0.0396    -2.48  2.11e- 2 TRUE 

The quadratic trend term is not statistically significant. Nevertheless, we will still fit a reduced model with a quadratic trend but we will omit the non-signficant seasonal Fourier term, cos1.

# Quadratic trend with standardized time variable

apple_reduced_quad_lm1 <- apple_ts |>
  model(apple_reduced_quad1 = TSLM(log(revenue_billions) ~ std_t + I(std_t^2) +
    sin1 + cos2 )) # Note sin2 is omitted
apple_reduced_quad_lm1 |>
  tidy() |>
  mutate(sig = p.value < 0.05)

# A tibble: 5 × 7
  .model              term        estimate std.error statistic  p.value sig  
  <chr>               <chr>          <dbl>     <dbl>     <dbl>    <dbl> <lgl>
1 apple_reduced_quad1 (Intercept)   1.83      0.0614    29.9   1.71e-20 TRUE 
2 apple_reduced_quad1 std_t         1.16      0.0415    28.0   7.94e-20 TRUE 
3 apple_reduced_quad1 I(std_t^2)   -0.0158    0.0478    -0.330 7.44e- 1 FALSE
4 apple_reduced_quad1 sin1          0.217     0.0573     3.80  8.79e- 4 TRUE 
5 apple_reduced_quad1 cos2         -0.0981    0.0410    -2.39  2.49e- 2 TRUE 

The quadratic trend term is still not statistically significant. We will fit a full model with a linear trend.

# Linear trend with standardized time variable

apple_full_linear_lm <- apple_ts |>
  model(apple_full_linear = TSLM(log(revenue_billions) ~ std_t + 
    sin1 + cos1 + cos2 )) # Note sin2 is omitted
apple_full_linear_lm |>
  tidy() |>
  mutate(sig = p.value < 0.05)

# A tibble: 5 × 7
  .model            term        estimate std.error statistic  p.value sig  
  <chr>             <chr>          <dbl>     <dbl>     <dbl>    <dbl> <lgl>
1 apple_full_linear (Intercept)   1.82      0.0388     46.9  4.04e-25 TRUE 
2 apple_full_linear std_t         1.16      0.0396     29.2  2.81e-20 TRUE 
3 apple_full_linear sin1          0.215     0.0540      3.99 5.42e- 4 TRUE 
4 apple_full_linear cos1          0.0934    0.0559      1.67 1.08e- 1 FALSE
5 apple_full_linear cos2         -0.0972    0.0388     -2.50 1.95e- 2 TRUE 

The coefficient on the cos1 seasonal Fourier term is not statistically significant. We now fit a reduced model that only contains the significant terms from the full model with a linear trend.

# Linear trend with standardized time variable

apple_reduced_linear_lm1 <- apple_ts |>
  model(apple_reduced_linear1 = TSLM(log(revenue_billions) ~ std_t + 
    sin1 + cos2 )) # Note sin2 is omitted
apple_reduced_linear_lm1 |>
  tidy() |>
  mutate(sig = p.value < 0.05)

# A tibble: 4 × 7
  .model                term        estimate std.error statistic  p.value sig  
  <chr>                 <chr>          <dbl>     <dbl>     <dbl>    <dbl> <lgl>
1 apple_reduced_linear1 (Intercept)   1.82      0.0402     45.3  1.60e-25 TRUE 
2 apple_reduced_linear1 std_t         1.16      0.0408     28.5  1.42e-20 TRUE 
3 apple_reduced_linear1 sin1          0.215     0.0559      3.85 7.21e- 4 TRUE 
4 apple_reduced_linear1 cos2         -0.0972    0.0402     -2.42 2.32e- 2 TRUE 

All the terms are statistically significant in this model. We now compare the models we have fitted using the AIC, AICc, and BIC criterion.

model_combined <- apple_ts |>
  model(
    apple_full_cubic = TSLM(log(revenue_billions) ~ std_t + I(std_t^2) + I(std_t^3) +
                    sin1 + cos1 + cos2),
    apple_full_quad = TSLM(log(revenue_billions) ~ std_t + I(std_t^2) + 
                    sin1 + cos1 + cos2),
    apple_reduced_quad1 = TSLM(log(revenue_billions) ~ std_t + I(std_t^2) + 
                    sin1 + cos2),
    apple_full_linear = TSLM(log(revenue_billions) ~ std_t +
                    sin1 + cos1 + cos2),
    apple_reduced_linear1 = TSLM(log(revenue_billions) ~ std_t + 
                    sin1 + cos2 )
  )

glance(model_combined) |>
  select(.model, AIC, AICc, BIC)

Table 2: Comparison of the AIC, AICc, and BIC values for the models fitted to the logarithm of the simulated time series.

Model	AIC	AICc	BIC
apple_full_cubic	-84	-76.8	-73.1
apple_full_quad	-82.5	-77.2	-73
apple_reduced_quad1	-81.3	-77.5	-73.1
apple_full_linear	**-84.4**	-80.6	-76.2
apple_reduced_linear1	-83.2	**-80.6**	**-76.4**

We will apply the apple_reduced_linear1 model.

Using the Residuals to Determine the Appropriate Correction

The residuals of this model are illustrated in Figure 5.

apple_resid_df <- model_combined |> 
  dplyr::select(apple_reduced_linear1) |>
  residuals() |> 
  as_tibble() |> 
  dplyr::select(.resid) |>
  rename(x = .resid) 
  
apple_resid_df |>
  mutate(density = dnorm(x, mean(apple_resid_df$x), sd(apple_resid_df$x))) |>
  ggplot(aes(x = x)) +
    geom_histogram(aes(y = after_stat(density)),
        color = "white", fill = "#56B4E9", binwidth = 0.05) +
    geom_line(aes(x = x, y = density)) +
    theme_bw() +
    labs(
      x = "Values",
      y = "Frequency",
      title = "Histogram of Residuals from the Reduced Model with a Linear Trend"
    ) +
    theme(
      plot.title = element_text(hjust = 0.5)
    )

Figure 5: Histogram of the residuals from the reduced model with a linear trend component

Using the command skewness(apple_resid_df$x), we compute the skewness of these residuals as: -0.799. This number is not close to zero (it is between -1 and -0.5) indicating moderate negative skewness. We would therefore apply the empirical correction to our mean forecast values.

Applying the Correction Factor

Table 3 summarizes some of the corrected mean values. Note that in this particular case, the corrected values are very close to the uncorrected values.

apple_model_values <- model_combined |> 
  dplyr::select(apple_reduced_linear1) |>
  glance()

apple_model_check <- apple_model_values |>
  mutate(
    sigma = sqrt(sigma2),
    lognorm_cf = exp((1/2) * sigma2),
    empirical_cf = apple_reduced_linear_lm1 |>
      residuals() |>
      pull(.resid) |>
      exp() |>
      mean()) |>
  select(.model, r_squared, sigma2, sigma, lognorm_cf, empirical_cf)

apple_pred <- model_combined |> 
  dplyr::select(apple_reduced_linear1) |>
  forecast(new_data = apple_ts) |>
  mutate(.mean_correction = .mean * apple_model_check$empirical_cf) |>
  select(t, revenue_billions, .mean, .mean_correction)

apple_model_check

# A tibble: 1 × 6
  .model                r_squared sigma2 sigma lognorm_cf empirical_cf
  <chr>                     <dbl>  <dbl> <dbl>      <dbl>        <dbl>
1 apple_reduced_linear1     0.971 0.0466 0.216       1.02         1.02

apple_pred <- model_combined |> 
  dplyr::select(apple_reduced_linear1) |>
  forecast(new_data = apple_ts) |>
  mutate(.mean_correction = .mean * apple_model_check$empirical_cf) |>
  select(t, revenue_billions, .mean, .mean_correction)

Table 3: Fitted values for the model representing Apple’s quarterly revenue

t	revenue_billions	.mean	.mean_correction	dates
1	t(N(0.22, 0.057))	1.284	1.309	2005-01-01
2	t(N(-0.051, 0.054))	0.975	0.994	2005-04-01
3	t(N(0.064, 0.057))	1.096	1.118	2005-07-01
4	t(N(0.22, 0.053))	1.281	1.306	2005-10-01
⋮	⋮	⋮	⋮	⋮
28	t(N(3.5, 0.054))	33.885	34.542	2011-10-01
29	t(N(4, 0.057))	58.588	59.725	2012-01-01

Plotting the Fitted Values

These fitted values are illustrated in Figure 6.

apple_ts |>
  autoplot(.vars = revenue_billions) +
  geom_line(data = apple_pred, aes(x = dates, y = .mean_correction), color = "#56B4E9") +
    labs(
      x = "Quarter",
      y = "Revenue (Billions USD)",
      title = "Apple Revenue in Billions of U.S. Dollars"
    ) +
    theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5))

Figure 6: Apple Inc.’s quarterly revenue in billions of U.S. dollars through first quarter of 2012 (in black) and the fitted regression model (in blue)

This time series was used as an example. We are obviously not interested in forecasting future values using this model. However, this is an excellent example of real-world exponential growth in a time series with a seasonal component. Limiting factors prevent exponential growth from being sustainable in the long run. After 2012, the Apple quarterly revenues follow a different, but very impressive, model. This is illustrated in Figure 7.

apple_raw |>
  dplyr::select(dates, revenue_billions) |>
  as_tsibble(index = dates) |>
  autoplot(.vars = revenue_billions) +
  geom_line(
    data = apple_raw |> filter(dates >= my("Jan 2012")), 
    aes(x = dates, y = revenue_billions), 
    color = "#D55E00"
  ) +
  labs(
    x = "Quarter",
    y = "Revenue (Billions USD)",
    title = "Apple Revenue (in Billions USD)"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

Figure 7: Apple Inc.’s quarterly revenue in billions of U.S. dollars; values beginning with the first quarter of 2012 are shown in orange

Choose One of the Following Small-Group Activities (25 min)

Small-Group Activity: Retail Sales (All Other General Merchandise Stores)

The code below downloads and gives the time plot for the total monthly sales in the United States for retail stores with the NAICS category 45299, “All Other General Merchandise Stores.” The time plot is given in Figure Figure 8.

# Read in retail sales data for "all other general merchandise stores"
retail_ts <- rio::import("https://byuistats.github.io/timeseries/data/retail_by_business_type.parquet") |>
  filter(naics == 45299) |>
  filter(as_date(month) >= my("Jan 1998")) |>
  mutate(t = 1:n()) |>
  mutate(std_t = (t - mean(t)) / sd(t)) |>
  mutate(
    cos1 = cos(2 * pi * 1 * t / 12),
    cos2 = cos(2 * pi * 2 * t / 12),
    cos3 = cos(2 * pi * 3 * t / 12),
    cos4 = cos(2 * pi * 4 * t / 12),
    cos5 = cos(2 * pi * 5 * t / 12),
    cos6 = cos(2 * pi * 6 * t / 12),
    sin1 = sin(2 * pi * 1 * t / 12),
    sin2 = sin(2 * pi * 2 * t / 12),
    sin3 = sin(2 * pi * 3 * t / 12),
    sin4 = sin(2 * pi * 4 * t / 12),
    sin5 = sin(2 * pi * 5 * t / 12)
  ) |>
  as_tsibble(index = month)

retail_ts |>
  autoplot(.vars = sales_millions) +
  stat_smooth(method = "lm", 
              formula = y ~ x, 
              geom = "smooth",
              se = FALSE,
              color = "#E69F00",
              linetype = "dotted") +
    labs(
      x = "Month",
      y = "Sales (Millions of U.S. Dollars)",
      title = paste0(retail_ts$business[1], " (", retail_ts$naics[1], ")")
    ) +
    theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5))

Figure 8: Time plot of the total monthly retail sales for all other general merchandise stores (45299)

Check Your Understanding

Use the retail sales data to do the following.

• Select an appropriate fitted model using the AIC, AICc, or BIC critera.
• Use the residuals to determine the appropriate correction for the data.
• Forecast the data for the next 5 years.
• Apply the appropriate correction to the forecasted values.
• Plot the fitted (forecasted) values along with the time series.

Small-Group Activity: Industrial Electricity Consumption in Texas

These data represent the amount of electricity used each month for industrial applications in Texas.

elec_ts <- rio::import("https://byuistats.github.io/timeseries/data/electricity_tx.csv") |>
  dplyr::select(-comments) |>
  mutate(month = my(month)) |>
  mutate(
    t = 1:n(),
    std_t = (t - mean(t)) / sd(t)
  ) |>
  mutate(
    cos1 = cos(2 * pi * 1 * t / 12),
    cos2 = cos(2 * pi * 2 * t / 12),
    cos3 = cos(2 * pi * 3 * t / 12),
    cos4 = cos(2 * pi * 4 * t / 12),
    cos5 = cos(2 * pi * 5 * t / 12),
    cos6 = cos(2 * pi * 6 * t / 12),
    sin1 = sin(2 * pi * 1 * t / 12),
    sin2 = sin(2 * pi * 2 * t / 12),
    sin3 = sin(2 * pi * 3 * t / 12),
    sin4 = sin(2 * pi * 4 * t / 12),
    sin5 = sin(2 * pi * 5 * t / 12)
  ) |>
  as_tsibble(index = month)

elec_plot_raw <- elec_ts |>
    autoplot(.vars = megawatthours) +
    labs(
      x = "Month",
      y = "Megawatt-hours",
      title = "Texas' Industrial Electricity Use"
    ) +
    theme_minimal() +
    theme(
      plot.title = element_text(hjust = 0.5)
    )

elec_plot_log <- elec_ts |>
    autoplot(.vars = log(megawatthours)) +
    labs(
      x = "Month",
      y = "log(Megwatt-hours)",
      title = "Log of Texas' Industrial Electricity Use"
    ) +
    theme_minimal() +
    theme(
      plot.title = element_text(hjust = 0.5)
    )

elec_plot_raw | elec_plot_log

Check Your Understanding

Use the Texas industrial electricity consumption data to do the following.

• Select an appropriate fitted model using the AIC, AICc, or BIC critera.
• Use the residuals to determine the appropriate correction for the data.
• Forecast the data for the next 5 years.
• Apply the appropriate correction to the forecasted values.
• Plot the fitted (forecasted) values along with the time series.

Homework Preview (5 min)

• Review upcoming homework assignment
• Clarify questions

Download Homework

homework_5_5.qmd

Small-Group Activity: Retail Sales

Check Your Understanding

Use the retail sales data to do the following.

• Select an appropriate fitted model using the AIC, AICc, or BIC critera.
• Use the residuals to determine the appropriate correction for the data.
• Forecast the data for the next 5 years.
• Apply the appropriate correction to the forecasted values.
• Plot the fitted (forecasted) values along with the time series.

retail_full_quad_lm <- retail_ts |>
  model(retail_full_quad = TSLM(log(sales_millions) ~ std_t + I(std_t^2) + 
    sin1 + cos1 + sin2 + cos2 + sin3 + cos3 
    + sin4 + cos4 + sin5 + cos5 + cos6 )) # Note sin6 is omitted
retail_full_quad_lm |>
  tidy() |>
  mutate(sig = p.value < 0.05) 

# A tibble: 14 × 7
   .model           term        estimate std.error statistic   p.value sig  
   <chr>            <chr>          <dbl>     <dbl>     <dbl>     <dbl> <lgl>
 1 retail_full_quad (Intercept)   8.37     0.00352   2378.   0         TRUE 
 2 retail_full_quad std_t         0.369    0.00235    157.   2.37e-279 TRUE 
 3 retail_full_quad I(std_t^2)    0.0130   0.00263      4.94 1.30e-  6 TRUE 
 4 retail_full_quad sin1         -0.0231   0.00332     -6.96 2.31e- 11 TRUE 
 5 retail_full_quad cos1          0.0353   0.00332     10.7  1.54e- 22 TRUE 
 6 retail_full_quad sin2         -0.0692   0.00332    -20.9  2.28e- 59 TRUE 
 7 retail_full_quad cos2          0.0785   0.00332     23.7  2.33e- 69 TRUE 
 8 retail_full_quad sin3         -0.0434   0.00332    -13.1  5.77e- 31 TRUE 
 9 retail_full_quad cos3          0.0699   0.00332     21.1  3.96e- 60 TRUE 
10 retail_full_quad sin4         -0.0288   0.00332     -8.69 2.83e- 16 TRUE 
11 retail_full_quad cos4          0.0627   0.00332     18.9  2.84e- 52 TRUE 
12 retail_full_quad sin5          0.0118   0.00332      3.57 4.25e-  4 TRUE 
13 retail_full_quad cos5          0.0636   0.00332     19.2  2.81e- 53 TRUE 
14 retail_full_quad cos6          0.0251   0.00235     10.7  1.04e- 22 TRUE 

retail_resid_df <- retail_full_quad_lm |> 
  residuals() |> 
  as_tibble() |> 
  dplyr::select(.resid) |>
  rename(x = .resid) 
  
retail_resid_df |>
  mutate(density = dnorm(x, mean(retail_resid_df$x), sd(retail_resid_df$x))) |>
  ggplot(aes(x = x)) +
    geom_histogram(aes(y = after_stat(density)),
        color = "white", fill = "#56B4E9", binwidth = 0.02) +
    geom_line(aes(x = x, y = density)) +
    theme_bw() +
    labs(
      x = "Values",
      y = "Frequency",
      title = "Histogram of Residuals from the Full Quadratic Model"
    ) +
    theme(
      plot.title = element_text(hjust = 0.5)
    )

skewness(retail_resid_df$x)

[1] 0.2205956

There is little skewness. We will use the log-normal correction factor.

Small-Group Activity: Texas Industrial Electricity Usage

Check Your Understanding

Use the Texas industrial electricity consumption data to do the following.

• Select an appropriate fitted model using the AIC, AICc, or BIC critera.
• Use the residuals to determine the appropriate correction for the data.
• Forecast the data for the next 5 years.
• Apply the appropriate correction to the forecasted values.
• Plot the fitted (forecasted) values along with the time series.

# Cubic model with standardized time variable

elec_full_cubic_lm <- elec_ts |>
  model(elec_full_cubic = TSLM(log(megawatthours) ~ std_t + I(std_t^2) + I(std_t^3) +
    sin1 + cos1 + sin2 + cos2 + sin3 + cos3 + sin4 + cos4 + sin5 + cos5 + cos6)) # Note sin6 is omitted
elec_full_cubic_lm |>
  tidy() |>
  mutate(sig = p.value < 0.05)

# A tibble: 15 × 7
   .model          term          estimate std.error statistic  p.value sig  
   <chr>           <chr>            <dbl>     <dbl>     <dbl>    <dbl> <lgl>
 1 elec_full_cubic (Intercept) 16.1         0.00684 2349.     0        TRUE 
 2 elec_full_cubic std_t        0.124       0.0115    10.8    1.19e-20 TRUE 
 3 elec_full_cubic I(std_t^2)   0.0313      0.00513    6.10   8.50e- 9 TRUE 
 4 elec_full_cubic I(std_t^3)   0.0105      0.00589    1.79   7.55e- 2 FALSE
 5 elec_full_cubic sin1        -0.0541      0.00648   -8.36   3.68e-14 TRUE 
 6 elec_full_cubic cos1        -0.0468      0.00645   -7.25   1.89e-11 TRUE 
 7 elec_full_cubic sin2        -0.0000864   0.00646   -0.0134 9.89e- 1 FALSE
 8 elec_full_cubic cos2        -0.000831    0.00645   -0.129  8.98e- 1 FALSE
 9 elec_full_cubic sin3         0.00103     0.00645    0.160  8.73e- 1 FALSE
10 elec_full_cubic cos3         0.0131      0.00645    2.03   4.46e- 2 TRUE 
11 elec_full_cubic sin4         0.0121      0.00645    1.87   6.29e- 2 FALSE
12 elec_full_cubic cos4         0.00272     0.00645    0.421  6.74e- 1 FALSE
13 elec_full_cubic sin5         0.0220      0.00645    3.41   8.36e- 4 TRUE 
14 elec_full_cubic cos5        -0.000505    0.00645   -0.0783 9.38e- 1 FALSE
15 elec_full_cubic cos6        -0.00172     0.00456   -0.378  7.06e- 1 FALSE

# Quadratic model with standardized time variable

elec_full_quadratic_lm <- elec_ts |>
  model(elec_full_cubic = TSLM(log(megawatthours) ~ std_t + I(std_t^2) + 
    sin1 + cos1 + sin2 + cos2 + sin3 + cos3 + sin4 + cos4 + sin5 + cos5 + cos6)) # Note sin6 is omitted
elec_full_quadratic_lm |>
  tidy() |>
  mutate(sig = p.value < 0.05)

# A tibble: 14 × 7
   .model          term         estimate std.error statistic  p.value sig  
   <chr>           <chr>           <dbl>     <dbl>     <dbl>    <dbl> <lgl>
 1 elec_full_cubic (Intercept) 16.1        0.00689 2333.     0        TRUE 
 2 elec_full_cubic std_t        0.143      0.00462   31.0    5.06e-68 TRUE 
 3 elec_full_cubic I(std_t^2)   0.0313     0.00517    6.05   1.04e- 8 TRUE 
 4 elec_full_cubic sin1        -0.0551     0.00650   -8.47   1.79e-14 TRUE 
 5 elec_full_cubic cos1        -0.0465     0.00650   -7.16   3.03e-11 TRUE 
 6 elec_full_cubic sin2        -0.000536   0.00650   -0.0825 9.34e- 1 FALSE
 7 elec_full_cubic cos2        -0.000571   0.00650   -0.0880 9.30e- 1 FALSE
 8 elec_full_cubic sin3         0.000775   0.00650    0.119  9.05e- 1 FALSE
 9 elec_full_cubic cos3         0.0133     0.00650    2.05   4.19e- 2 TRUE 
10 elec_full_cubic sin4         0.0119     0.00649    1.84   6.80e- 2 FALSE
11 elec_full_cubic cos4         0.00298    0.00650    0.458  6.48e- 1 FALSE
12 elec_full_cubic sin5         0.0219     0.00649    3.37   9.39e- 4 TRUE 
13 elec_full_cubic cos5        -0.000245   0.00650   -0.0378 9.70e- 1 FALSE
14 elec_full_cubic cos6        -0.00159    0.00459   -0.347  7.29e- 1 FALSE

elec_resid_df <- elec_full_quadratic_lm |> 
  residuals() |> 
  as_tibble() |> 
  dplyr::select(.resid) |>
  rename(x = .resid) 
  
elec_resid_df |>
  mutate(density = dnorm(x, mean(elec_resid_df$x), sd(elec_resid_df$x))) |>
  ggplot(aes(x = x)) +
    geom_histogram(aes(y = after_stat(density)),
        color = "white", fill = "#56B4E9", binwidth = 0.02) +
    geom_line(aes(x = x, y = density)) +
    theme_bw() +
    labs(
      x = "Values",
      y = "Frequency",
      title = "Histogram of Residuals from the Full Quadratic Model"
    ) +
    theme(
      plot.title = element_text(hjust = 0.5)
    )

skewness(elec_resid_df$x)

[1] -0.9556189

There is moderate negative skewness. We will use the emperical correction factor.