Half-hourly electricity demand for Victoria, Australia is contained in
vic_elec. Extract the January 2014 electricity demand, and aggregate this data to daily with daily total demands and maximum temperatures.jan14_vic_elec <- vic_elec |> filter(yearmonth(Time) == yearmonth("2014 Jan")) |> index_by(Date = as_date(Time)) |> summarise(Demand = sum(Demand), Temperature = max(Temperature))
- Plot the data and find the regression model for Demand with temperature as a predictor variable. Why is there a positive relationship?
jan14_vic_elec <- vic_elec |>
filter(yearmonth(Time) == yearmonth("2014 Jan")) |>
index_by(Date = as_date(Time)) |>
summarise(Demand = sum(Demand), Temperature = max(Temperature))
jan14_vic_elec |>
pivot_longer(2:3, names_to = "key", values_to = "value") |>
autoplot(.vars = value) +
facet_grid(vars(key), scales = "free_y")
jan14_vic_elec |>
ggplot(aes(x = Temperature, y = Demand)) +
geom_point()
fit <- jan14_vic_elec |>
model(TSLM(Demand ~ Temperature))
fit |> report()Series: Demand
Model: TSLM
Residuals:
Min 1Q Median 3Q Max
-49978.2 -10218.9 -121.3 18533.2 35440.6
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 59083.9 17424.8 3.391 0.00203 **
Temperature 6154.3 601.3 10.235 3.89e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 24540 on 29 degrees of freedom
Multiple R-squared: 0.7832, Adjusted R-squared: 0.7757
F-statistic: 104.7 on 1 and 29 DF, p-value: 3.8897e-11There is a clear positive relationship: each extra degree of maximum temperature adds about 6,150 MWh to daily demand. January is summer in Victoria, so hotter days lead to more electricity demand (presumably for air conditioning).
- Produce a residual plot. Is the model adequate? Are there any outliers or influential observations?
fit |> gg_tsresiduals()
The residuals tend to rise over the month, suggesting that the model is lacking a trend (although looking at a larger window of the data, this is clearly not a real pattern). The largest (negative) residuals are on public holidays and weekends (1 January, Sunday 5 January, and 26–27 January around Australia Day) when demand is lower than the temperature alone would suggest. These are outliers, although none appears particularly influential. The hottest days (above 40^\circ\text{C}) have the most leverage, but they lie close to the fitted line. The lag 1 autocorrelation of the residuals (0.34) is just inside the critical bounds.
Use the model to forecast the electricity demand that you would expect for the next day if the maximum temperature was 15^\circ\text{C} and compare it with the forecast if the with maximum temperature was 35^\circ\text{C}. Do you believe these forecasts? The following R code will get you started:
jan14_vic_elec |> model(TSLM(Demand ~ Temperature)) |> forecast(new_data(jan14_vic_elec, 1) |> mutate(Temperature = 15)) |> autoplot(jan14_vic_elec)
next_day <- scenarios(
`Cold day` = new_data(jan14_vic_elec, 1) |> mutate(Temperature = 15),
`Hot day` = new_data(jan14_vic_elec, 1) |> mutate(Temperature = 35)
)
fc <- fit |>
forecast(new_data = next_day)
autoplot(jan14_vic_elec, Demand) +
autolayer(fc[1, ]) +
autolayer(fc[2, ])
The forecast is 151,398 MWh for a 15^\circ\text{C} day and 274,484 MWh for a 35^\circ\text{C} day. The forecast for 35^\circ\text{C} seems reasonable, as this temperature is well within the range of the January data. The forecast for 15^\circ\text{C} is not believable: the lowest daily maximum observed in January 2014 was 19.6^\circ\text{C} (January is a summer month in Victoria), so this is an extrapolation beyond the data. As part (e) shows, demand actually rises on cold days (due to heating), so this forecast is likely to be much too low.
- Give prediction intervals for your forecasts.
fc |>
hilo() |>
select(-.model)# A tsibble: 2 x 7 [1D]
# Key: .scenario [2]
.scenario Date
<chr> <date>
1 Cold day 2014-02-01
2 Hot day 2014-02-01
# ℹ 5 more variables: Demand <dist>, .mean <dbl>, Temperature <dbl>,
# `80%` <hilo>, `95%` <hilo>The intervals are wide (about \pm 50{,}000 MWh for the 95% intervals), and they assume the linear model is correct. That is doubtful for the 15^\circ\text{C} forecast, which is outside the range of the data.
- Plot Demand vs Temperature for all of the available data in
vic_elecaggregated to daily total demand and maximum temperature. What does this say about your model?
vic_elec |>
index_by(Date = as_date(Time)) |>
summarise(Demand = sum(Demand), Temperature = max(Temperature)) |>
ggplot(aes(x = Temperature, y = Demand)) +
geom_point()
The “U” shaped scatterplot suggests that there is a non-linear relationship between electricity demand and daily maximum temperature. The model assumption that electricity demand is linearly predicted by temperature is incorrect.
Data set
olympic_runningcontains the winning times (in seconds) in each Olympic Games sprint, middle-distance and long-distance track events from 1896 to 2016.
- Plot the winning time against the year for each event. Describe the main features of the plot.
olympic_running |>
ggplot(aes(x = Year, y = Time, colour = Sex)) +
geom_line() +
geom_point(size = 1) +
facet_wrap(~Length, scales = "free_y", nrow = 2) +
theme_minimal() +
scale_color_brewer(palette = "Dark2") +
theme(legend.position = "bottom", legend.title = element_blank()) +
labs(y = "Running time (seconds)")
The running times are generally decreasing as time progresses (although the rate of this decline is slowing down in recent Olympics). There are some missing values in the data corresponding to the World Wars (in which the Olympic Games were not held), and the women’s events started much later than the men’s.
- Fit a regression line to the data for each event. Obviously the winning times have been decreasing, but at what average rate per year?
fit <- olympic_running |>
model(TSLM(Time ~ trend()))
tidy(fit)# A tibble: 28 × 8
Length Sex .model term estimate std.error statistic p.value
<int> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 100 men TSLM(Time ~ trend()) (Int… 11.1 0.0909 123. 1.86e-37
2 100 men TSLM(Time ~ trend()) tren… -0.0504 0.00479 -10.5 7.24e-11
3 100 women TSLM(Time ~ trend()) (Int… 12.4 0.163 76.0 4.58e-25
4 100 women TSLM(Time ~ trend()) tren… -0.0567 0.00749 -7.58 3.72e- 7
5 200 men TSLM(Time ~ trend()) (Int… 22.3 0.140 159. 4.06e-39
6 200 men TSLM(Time ~ trend()) tren… -0.0995 0.00725 -13.7 3.80e-13
7 200 women TSLM(Time ~ trend()) (Int… 25.5 0.525 48.6 8.34e-19
8 200 women TSLM(Time ~ trend()) tren… -0.135 0.0227 -5.92 2.17e- 5
9 400 men TSLM(Time ~ trend()) (Int… 50.3 0.445 113. 1.53e-36
10 400 men TSLM(Time ~ trend()) tren… -0.258 0.0235 -11.0 2.75e-11
# ℹ 18 more rowsThe Games are held every four years, so the trend() coefficient is the average change per Games, and we divide it by 4 to get the average change per year.
The men’s 100 running time has been decreasing by an average of 0.013 seconds each year.
The women’s 100 running time has been decreasing by an average of 0.014 seconds each year.
The men’s 200 running time has been decreasing by an average of 0.025 seconds each year.
The women’s 200 running time has been decreasing by an average of 0.034 seconds each year.
The men’s 400 running time has been decreasing by an average of 0.065 seconds each year.
The women’s 400 running time has been decreasing by an average of 0.04 seconds each year.
The men’s 800 running time has been decreasing by an average of 0.152 seconds each year.
The women’s 800 running time has been decreasing by an average of 0.198 seconds each year.
The men’s 1500 running time has been decreasing by an average of 0.315 seconds each year.
The women’s 1500 running time has been increasing by an average of 0.147 seconds each year.
The men’s 5000 running time has been decreasing by an average of 1.03 seconds each year.
The women’s 5000 running time has been decreasing by an average of 0.303 seconds each year.
The men’s 10000 running time has been decreasing by an average of 2.666 seconds each year.
The women’s 10000 running time has been decreasing by an average of 3.496 seconds each year.
The slopes give the average change in winning time per year: for example, 0.013 seconds per year for the men’s 100m and 2.67 seconds per year for the men’s 10000m. All the winning times have been decreasing except the women’s 1500m, where the times since the event was introduced in 1972 have been getting slower on average (the fastest times were in the 1980s and early 1990s, and the 2012 and 2016 races were the slowest).
- Plot the residuals against the year. What does this indicate about the suitability of the fitted lines?
augment(fit) |>
ggplot(aes(x = Year, y = .innov, colour = Sex)) +
geom_line() +
geom_point(size = 1) +
facet_wrap(~Length, scales = "free_y", nrow = 2) +
theme_minimal() +
scale_color_brewer(palette = "Dark2") +
theme(legend.position = "bottom", legend.title = element_blank())
It doesn’t seem that the linear trend is appropriate for this data. For many of the men’s events, the residuals show curvature: they are mostly positive at the start and end of the series, and negative in the middle. This reflects a change in slope: the winning times fell quickly in the early decades, but have improved much more slowly in recent decades. There are also some outliers, most obviously the 1896 times for the men’s 400m, 800m and 1500m, which were much slower than the fitted lines.
- Predict the winning time for each race in the 2020 Olympics. Give a prediction interval for your forecasts. What assumptions have you made in these calculations?
fit |>
forecast(h = 1) |>
mutate(PI = hilo(Time, 95)) |>
select(-.model)# A fable: 14 x 6 [4Y]
# Key: Length, Sex [14]
Length Sex Year
<int> <chr> <dbl>
1 100 men 2020
2 100 women 2020
3 200 men 2020
4 200 women 2020
5 400 men 2020
6 400 women 2020
7 800 men 2020
8 800 women 2020
9 1500 men 2020
10 1500 women 2020
11 5000 men 2020
12 5000 women 2020
13 10000 men 2020
14 10000 women 2020
# ℹ 3 more variables: Time <dist>, .mean <dbl>, PI <hilo>Using a linear trend we assume that winning times will decrease in a linear fashion which is unrealistic for running times. As we saw from the residual plots above there are mainly large positive residuals for the last few years, indicating that the decreases in the winning times are not linear. Because the rate of improvement has slowed, assuming a constant rate means the forecasts are probably too low (i.e., too optimistic). We also assume that the residuals are normally distributed and uncorrelated with constant variance.
The 2020 Olympics were held in Tokyo in 2021, so we can compare the forecasts with the actual winning times.
tokyo <- tribble(
~Length, ~Sex, ~Actual,
100, "men", 9.80, 100, "women", 10.61,
200, "men", 19.62, 200, "women", 21.53,
400, "men", 43.85, 400, "women", 48.36,
800, "men", 105.06, 800, "women", 115.21,
1500, "men", 208.32, 1500, "women", 233.11,
5000, "men", 778.15, 5000, "women", 876.79,
10000, "men", 1663.22, 10000, "women", 1795.32
)
fit |>
forecast(h = 1) |>
hilo(95) |>
as_tibble() |>
select(Length, Sex, .mean, `95%`) |>
left_join(tokyo, by = c("Length", "Sex"))# A tibble: 14 × 5
Length Sex .mean `95%` Actual
<dbl> <chr> <dbl> <hilo> <dbl>
1 100 men 9.53 [ 9.049446, 10.01957]95 9.8
2 100 women 10.5 [ 10.056992, 11.01183]95 10.6
3 200 men 19.1 [ 18.423087, 19.81179]95 19.6
4 200 women 21.2 [ 20.106551, 22.29293]95 21.5
5 400 men 42.0 [ 39.668596, 44.41603]95 43.8
6 400 women 48.4 [ 46.123295, 50.74923]95 48.4
7 800 men 99.2 [ 92.137378, 106.33802]95 105.
8 800 women 111. [ 104.105124, 118.86174]95 115.
9 1500 men 207. [ 190.091940, 223.90462]95 208.
10 1500 women 245. [ 233.786710, 257.12632]95 233.
11 5000 men 773. [ 739.735488, 805.91363]95 778.
12 5000 women 892. [ 814.664722, 969.56461]95 877.
13 10000 men 1580. [1519.300480, 1641.40196]95 1663.
14 10000 women 1763. [1717.895285, 1808.13543]95 1795. The actual winning time was slower than the point forecast in 11 of the 14 events, consistent with the forecasts being too optimistic. For example, the men’s 400m was won in 43.85 seconds, well above the point forecast of 42.04 seconds, although inside the 95% prediction interval. Two of the actual times lay outside the 95% prediction intervals: the men’s 10000m was much slower than predicted, while the women’s 1500m was faster.
An elasticity coefficient is the ratio of the percentage change in the forecast variable (y) to the percentage change in the predictor variable (x). Mathematically, the elasticity is defined as (dy/dx) \times (x/y). Consider the log-log model, \log y = \beta_0 + \beta_1 \log x + \varepsilon. Express y as a function of x and show that the coefficient \beta_1 is the elasticity coefficient.
Exponentiating both sides gives y = e^{\beta_0 + \varepsilon} x^{\beta_1}. Then dy/dx = \beta_1 e^{\beta_0 + \varepsilon} x^{\beta_1 - 1} = \beta_1 y/x, so (dy/dx) \times (x/y) = \beta_1.
Equivalently, differentiating both sides of the log-log model with respect to x (holding \varepsilon fixed) gives \frac{1}{y}\frac{dy}{dx} = \frac{\beta_1}{x}, so \beta_1 = \frac{dy}{dx}\times\frac{x}{y}.
The data set
souvenirsconcerns the monthly sales figures of a shop which opened in January 1987 and sells gifts, souvenirs, and novelties. The shop is situated on the wharf at a beach resort town in Queensland, Australia. The sales volume varies with the seasonal population of tourists. There is a large influx of visitors to the town at Christmas and for the local surfing festival, held every March since 1988. Over time, the shop has expanded its premises, range of products, and staff.
- Produce a time plot of the data and describe the patterns in the graph. Identify any unusual or unexpected fluctuations in the time series.
souvenirs |> autoplot(Sales)
Features of the data:
- Explain why it is necessary to take logarithms of these data before fitting a model.
The last feature above suggests taking logs to make the pattern (and variance) more stable.
# Taking logarithms of the data
souvenirs |> autoplot(log(Sales))
After taking logs, the trend looks more linear and the seasonal variation is roughly constant.
- Fit a regression model to the logarithms of these sales data with a linear trend, seasonal dummies and a “surfing festival” dummy variable.
fit <- souvenirs |>
mutate(festival = month(Month) == 3 & year(Month) != 1987) |>
model(reg = TSLM(log(Sales) ~ trend() + season() + festival))
souvenirs |>
autoplot(Sales, col = "gray") +
geom_line(data = augment(fit), aes(y = .fitted), col = "blue")
- Plot the residuals against time and against the fitted values. Do these plots reveal any problems with the model?
fit |> gg_tsresiduals()
The residuals are serially correlated. This is visible in both the time plot and the ACF. The residuals reveal nonlinearity in the trend.
augment(fit) |>
ggplot(aes(x = .fitted, y = .innov)) +
geom_point() +
scale_x_log10()
The plot of residuals against fitted values looks fine – no notable patterns emerge. We plot the fitted values on a log scale because we took logs in the model.
- Do boxplots of the residuals for each month. Does this reveal any problems with the model?
augment(fit) |>
mutate(month = month(Month, label = TRUE)) |>
ggplot(aes(x = month, y = .innov)) +
geom_boxplot()
The boxplots show differences in variation across the months revealing some potential heteroscedasticity.
- What do the values of the coefficients tell you about each variable?
tidy(fit) |> mutate(pceffect = (exp(estimate) - 1) * 100)# A tibble: 14 × 7
.model term estimate std.error statistic p.value pceffect
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 reg (Intercept) 7.62 0.0742 103. 4.67e-78 203688.
2 reg trend() 0.0220 0.000827 26.6 2.32e-38 2.23
3 reg season()year2 0.251 0.0957 2.63 1.06e- 2 28.6
4 reg season()year3 0.266 0.193 1.38 1.73e- 1 30.5
5 reg season()year4 0.384 0.0957 4.01 1.48e- 4 46.8
6 reg season()year5 0.409 0.0957 4.28 5.88e- 5 50.6
7 reg season()year6 0.449 0.0958 4.69 1.33e- 5 56.6
8 reg season()year7 0.610 0.0958 6.37 1.71e- 8 84.1
9 reg season()year8 0.588 0.0959 6.13 4.53e- 8 80.0
10 reg season()year9 0.669 0.0959 6.98 1.36e- 9 95.3
11 reg season()year10 0.747 0.0960 7.79 4.48e-11 111.
12 reg season()year11 1.21 0.0960 12.6 1.29e-19 234.
13 reg season()year12 1.96 0.0961 20.4 3.39e-31 612.
14 reg festivalTRUE 0.502 0.196 2.55 1.29e- 2 65.1 Because the model is for log sales, 100(e^{\beta} - 1) is the percentage effect of each variable on sales.
(Intercept) is not interpretable.trend() coefficient shows that with every month sales increase on average by 2.2%.season()year2 coefficient shows that February sales are greater than January on average by 28.6%, after allowing for the trend.season()year12 coefficient shows that December sales are greater than January on average by 611.5%, after allowing for the trend.festivalTRUE coefficient shows that for months that include the surfing festival, sales increase on average by 65.1% compared to months without the festival, after allowing for the trend and seasonality.All the coefficients are significant at the 5% level except season()year3: once the festival is allowed for, March sales are not significantly different from January sales. (As the festival has been held every March since 1988, season()year3 is estimated only from March 1987.)
- What does the Ljung-Box test tell you about your model?
augment(fit) |>
features(.innov, ljung_box, lag = 24)# A tibble: 1 × 3
.model lb_stat lb_pvalue
<chr> <dbl> <dbl>
1 reg 112. 2.15e-13The serial correlation in the residuals is significant (the p-value is tiny). So the prediction intervals will be too narrow, and there is information in the residuals that the model has not used.
- Regardless of your answers to the above questions, use your regression model to predict the monthly sales for 1994, 1995, and 1996. Produce prediction intervals for each of your forecasts.
The festival is held every March in the forecast period. The forecasts and prediction intervals are computed on the log scale, and automatically back-transformed.
future_souvenirs <- new_data(souvenirs, n = 36) |>
mutate(festival = month(Month) == 3)
fit |>
forecast(new_data = future_souvenirs) |>
autoplot(souvenirs)
- How could you improve these predictions by modifying the model?
The
us_gasolineseries consists of weekly data for supplies of US finished motor gasoline product, from 2 February 1991 to 20 January 2017. The units are in “million barrels per day”. Consider only the data to the end of 2004.
- Fit a harmonic regression with trend to the data. Experiment with changing the number Fourier terms. Plot the observed gasoline and fitted values and comment on what you see.
gas <- us_gasoline |> filter(year(Week) <= 2004)
gas |> autoplot(Barrels)
We fit models with K = 1, 2, \dots, 15 pairs of Fourier terms, and compute the AICc and CV statistics for each. For weekly data, fourier() uses a seasonal period of 52 weeks (rather than the average of 365.25/7 \approx 52.18 weeks in a year).
fit <- gas |>
model(
fourier1 = TSLM(Barrels ~ trend() + fourier(K = 1)),
fourier2 = TSLM(Barrels ~ trend() + fourier(K = 2)),
fourier3 = TSLM(Barrels ~ trend() + fourier(K = 3)),
fourier4 = TSLM(Barrels ~ trend() + fourier(K = 4)),
fourier5 = TSLM(Barrels ~ trend() + fourier(K = 5)),
fourier6 = TSLM(Barrels ~ trend() + fourier(K = 6)),
fourier7 = TSLM(Barrels ~ trend() + fourier(K = 7)),
fourier8 = TSLM(Barrels ~ trend() + fourier(K = 8)),
fourier9 = TSLM(Barrels ~ trend() + fourier(K = 9)),
fourier10 = TSLM(Barrels ~ trend() + fourier(K = 10)),
fourier11 = TSLM(Barrels ~ trend() + fourier(K = 11)),
fourier12 = TSLM(Barrels ~ trend() + fourier(K = 12)),
fourier13 = TSLM(Barrels ~ trend() + fourier(K = 13)),
fourier14 = TSLM(Barrels ~ trend() + fourier(K = 14)),
fourier15 = TSLM(Barrels ~ trend() + fourier(K = 15))
)
glance(fit) |>
arrange(AICc) |>
select(.model, AICc, CV)# A tibble: 15 × 3
.model AICc CV
<chr> <dbl> <dbl>
1 fourier7 -1887. 0.0740
2 fourier11 -1885. 0.0742
3 fourier8 -1885. 0.0743
4 fourier12 -1884. 0.0742
5 fourier6 -1884. 0.0744
6 fourier9 -1882. 0.0745
7 fourier10 -1881. 0.0746
8 fourier13 -1881. 0.0745
9 fourier14 -1879. 0.0748
10 fourier15 -1876. 0.0750
11 fourier5 -1862. 0.0766
12 fourier4 -1856. 0.0773
13 fourier3 -1853. 0.0777
14 fourier2 -1814. 0.0819
15 fourier1 -1809. 0.0825# Best model has order 7
gas |>
autoplot(Barrels, col = "gray") +
geom_line(
data = augment(fit) |> filter(.model == "fourier7"),
aes(y = .fitted),
col = "blue"
)
Increasing K reduces the AICc and CV quickly up to K=7, after which they change very little, and both increase again for larger K. The fitted values from the best model (selected in part b) seem to have captured the annual seasonality quite well.
- Select the appropriate number of Fourier terms to include by minimising the AICc or CV value.
Both the AICc and the CV statistic are minimised with K = 7, as shown in the table above (so the model has 14 Fourier terms). The values for K between 6 and 12 are very close, so a model with K=6 would give very similar results.
- Plot the residuals of the final model using the
gg_tsresiduals()function and comment on these. Use a Ljung-Box test to check for residual autocorrelation.
fit |>
select(fourier7) |>
gg_tsresiduals()
There is some small remaining autocorrelation at lag 1. The residuals look roughly normal, with constant variance.
augment(fit) |>
filter(.model == "fourier7") |>
features(.innov, ljung_box, lag = 52)# A tibble: 1 × 3
.model lb_stat lb_pvalue
<chr> <dbl> <dbl>
1 fourier7 56.9 0.299The Ljung-Box test is not significant (p-value = 0.30). Even if it had given a significant result here, the correlations are very small, so it would not make much difference to the forecasts and prediction intervals.
- Generate forecasts for the next year of data and plot these along with the actual data for 2005. Comment on the forecasts.
fit |>
select("fourier7") |>
forecast(h = "1 year") |>
autoplot(filter_index(us_gasoline, "2000" ~ "2006"))
The forecasts fit pretty well for the first eight months, but not so well after that, when supplies fell sharply following Hurricane Katrina at the end of August 2005.
The annual population of Afghanistan is available in the
global_economydata set.
- Plot the data and comment on its features. Can you observe the effect of the Soviet-Afghan war?
global_economy |>
filter(Country == "Afghanistan") |>
autoplot(Population / 1e6) +
labs(y = "Population (millions)") +
geom_ribbon(aes(xmin = 1979.98, xmax = 1989.13), fill = "pink", alpha = 0.4) +
annotate(
"text",
x = 1984.5,
y = 10,
label = "Soviet-Afghan war",
col = "red",
size = 3
)
The population increases slowly from 1960 to 1980, then decreases during the Soviet-Afghan war (24 Dec 1979 – 15 Feb 1989), and has increased relatively rapidly since then. The last 30 years has shown an almost linear increase in population.
- Fit a linear trend model and compare this to a piecewise linear trend model with knots at 1980 and 1989.
fit <- global_economy |>
filter(Country == "Afghanistan") |>
model(
linear = TSLM(Population ~ trend()),
piecewise = TSLM(Population ~ trend(knots = c(1980, 1989)))
)
augment(fit) |>
autoplot(.fitted) +
geom_line(aes(y = Population), colour = "black")
The fitted values show that the piecewise linear model has tracked the data very closely, while the linear model is inaccurate.
- Generate forecasts from these two models for the five years after the end of the data, and comment on the results.
fc <- fit |> forecast(h = "5 years")
autoplot(fc) +
autolayer(
filter(global_economy |> filter(Country == "Afghanistan")),
Population
)
The linear model is clearly incorrect with prediction intervals too wide, and the point forecasts too low.
The piecewise linear model looks good, but the prediction intervals are probably too narrow. This model assumes that the trend since the last knot will continue unchanged, which is unrealistic.
(For advanced readers following on from Section 7.9.)
Using matrix notation it was shown that if \mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}, where \boldsymbol{\varepsilon} has mean \mathbf{0} and variance matrix \sigma^2\mathbf{I}, the estimated coefficients are given by \hat{\boldsymbol{\beta}} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{y} and a forecast is given by \hat{y} = \mathbf{x}^*\hat{\boldsymbol{\beta}} = \mathbf{x}^*(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{y} where \mathbf{x}^* is a row vector containing the values of the predictors for the forecast (in the same format as \mathbf{X}), and the forecast variance is given by \text{Var}(\hat{y}) = \hat{\sigma}^2 \left[ 1 + \mathbf{x}^*(\mathbf{X}'\mathbf{X})^{-1}(\mathbf{x}^*)' \right].
Consider the simple time trend model where y_t = \beta_0 + \beta_1t + \varepsilon_t. Using the following results, \sum^{T}_{t=1}{t} = \frac{1}{2}T(T + 1), \quad \sum^{T}_{t=1}{t^2} = \frac{1}{6}T(T + 1)(2T + 1) derive the following expressions:
- \displaystyle\mathbf{X}'\mathbf{X} = \frac{1}{6}\left[ \begin{array}{cc} 6T & 3T(T + 1) \\ 3T(T + 1) & T(T + 1)(2T + 1) \\ \end{array} \right]
- \displaystyle(\mathbf{X}'\mathbf{X})^{-1} = \frac{2}{T(T^2 - 1)}\left[ \begin{array}{cc} (T + 1)(2T + 1) &-3(T + 1) \\-3(T + 1) & 6 \\ \end{array} \right]
- \displaystyle\hat{\beta}_0 = \frac{2}{T(T - 1)}\left[ (2T + 1)\sum^T_{t=1}y_t - 3\sum^T_{t=1}ty_t \right]
\displaystyle\hat{\beta}_1 = \frac{6}{T(T^2 - 1)}\left[ 2\sum^T_{t=1}ty_t - (T + 1)\sum^T_{t=1}y_t \right]
- \displaystyle\text{Var}(\hat{y}_{t}) = \hat{\sigma}^2\left[ 1 + \frac{2}{T(T - 1)}\left( 1 - 4T - 6h + 6\frac{(T + h)^2}{T + 1} \right) \right]
For a time trend, \mathbf{X} = \begin{bmatrix} 1 & 1 \\ 1 & 2 \\ 1 & 3 \\ \vdots \\ 1 & T \end{bmatrix}
\begin{align*} \mathbf{X}' \mathbf{X} & = \begin{bmatrix} 1 & 1 & 1 & \dots & 1 \\ 1 & 2 & 3 & \dots & T \end{bmatrix} \begin{bmatrix} 1 & 1 \\ 1 & 2 \\ 1 & 3 \\ \vdots \\ 1 & T \end{bmatrix} \\ & = \begin{bmatrix} T & 1 + 2 + 3 + \cdots + T \\ 1 + 2 + 3 + \cdots + T & 1^2 + 2^2 + 3^2 + \cdots + T^2 \end{bmatrix} \\ & = \begin{bmatrix} T & \frac{1}{2}T(T + 1) \\ \frac{1}{2}T(T + 1) & \frac{1}{6}T(T + 1)(2T + 1) \end{bmatrix} \\ & = \frac{1}{6}\begin{bmatrix} 6T & 3T(T + 1) \\ 3T(T + 1) & T(T + 1)(2T + 1) \end{bmatrix} \end{align*}
\begin{align*} (\mathbf{X}' \mathbf{X} )^{-1} & = \frac{6}{(6T)[T(T + 1)(2T + 1)] - 9T^2(T + 1)^2} \begin{bmatrix} T(T + 1)(2T + 1) & -3T(T + 1) \\ -3T(T + 1) & 6T \end{bmatrix} \\ & = \frac{6}{6[T(T + 1)(2T + 1)] - 9T(T + 1)^2} \begin{bmatrix} (T + 1)(2T + 1) & -3(T + 1) \\ -3(T + 1) & 6 \end{bmatrix} \\ & = \frac{2}{2T(T + 1)(2T + 1) - 3T(T + 1)^2} \begin{bmatrix} (T + 1)(2T + 1) & -3(T + 1) \\ -3(T + 1) & 6 \end{bmatrix} \\ & = \frac{2}{T(T + 1)[2(2T + 1) - 3(T + 1)]} \begin{bmatrix} (T + 1)(2T + 1) & -3(T + 1) \\ -3(T + 1) & 6 \end{bmatrix} \\ & = \frac{2}{T(T + 1)(T - 1)} \begin{bmatrix} (T + 1)(2T + 1) & -3(T + 1) \\ -3(T + 1) & 6 \end{bmatrix} \\ & = \frac{2}{T(T^2 - 1)} \begin{bmatrix} (T + 1)(2T + 1) & -3(T + 1) \\ -3(T + 1) & 6 \end{bmatrix} \\ \end{align*}
\begin{align*} \mathbf{X}' \mathbf{Y} = \begin{bmatrix} 1 & 1 & 1 & \dots & 1 \\ 1 & 2 & 3 & \dots & T \end{bmatrix}\mathbf{Y} = \begin{bmatrix} \sum_{i=1}^T Y_i \\ \sum_{i=1}^T iY_i \end{bmatrix} \end{align*}
Now the first element of (\mathbf{X}' \mathbf{X} )^{-1} \mathbf{X}' \mathbf{Y} is \begin{align*} \hat{\beta}_0 & = \frac{2}{T(T^2 - 1)} \begin{bmatrix} (T + 1)(2T + 1) & -3(T + 1) \end{bmatrix} \begin{bmatrix} \sum_{i=1}^T Y_i \\ \sum_{i=1}^T iY_i \end{bmatrix} \\ & = \frac{2}{T(T^2 - 1)} \left[ (T + 1)(2T + 1)\sum_{i=1}^T Y_i - 3(T + 1) \sum_{i=1}^T iY_i \right] \\ & = \frac{2(T + 1)}{T(T^2 - 1)} \left[ (2T + 1)\sum_{i=1}^T Y_i - 3 \sum_{i=1}^T iY_i \right] \\ & = \frac{2}{T(T - 1)} \left[ (2T + 1)\sum_{i=1}^T Y_i - 3 \sum_{i=1}^T iY_i \right] . \end{align*}
The second element of (\mathbf{X}' \mathbf{X} )^{-1} \mathbf{X}' \mathbf{Y} is \begin{align*} \hat{\beta}_1 & = \frac{2}{T(T^2 - 1)} \begin{bmatrix} -3(T + 1) & 6 \end{bmatrix} \begin{bmatrix} \sum_{i=1}^T Y_i \\ \sum_{i=1}^T iY_i \end{bmatrix} \\ & = \frac{2}{T(T^2 - 1)} \left[ -3(T + 1) \sum_{i=1}^T Y_i + 6 \sum_{i=1}^T iY_i \right] \\ & = \frac{6}{T(T^2 - 1)} \left[ 2 \sum_{i=1}^T iY_i - (T + 1) \sum_{i=1}^T Y_i \right] \end{align*}
Now \mathbf{X}^* = [1 ~~~ T + h]. So \text{Var}(Y^*| \mathbf{Y},\mathbf{X},\mathbf{X}^*) \begin{align*} & = \hat{\sigma}^2\left\{ 1 + \mathbf{X}^*(\mathbf{X}'\mathbf{X})^{-1}(\mathbf{X}^*)' \right\} \\ & = \hat{\sigma}^2\left\{ 1 + \frac{2}{T(T^2 - 1)} \begin{bmatrix} 1 & T + h \end{bmatrix} \begin{bmatrix} (T + 1)(2T + 1) &-3(T + 1) \\-3(T + 1) & 6 \end{bmatrix} \begin{bmatrix} 1 \\ T + h \end{bmatrix} \right\} \\ & = \hat{\sigma}^2\left\{ 1 + \frac{2}{T(T^2 - 1)} \begin{bmatrix} (T + 1)(2T + 1) - 3(T + 1)(T + h) &-3(T + 1) + 6(T + h) \end{bmatrix} \begin{bmatrix} 1 \\ T + h \end{bmatrix} \right\} \\ & = \hat{\sigma}^2\left\{ 1 + \frac{2}{T(T^2 - 1)} \left[ (T + 1)(2T + 1) - 6(T + 1)(T + h) + 6(T + h)^2 \right] \right\} \\ & = \hat{\sigma}^2\left\{ 1 + \frac{2}{T(T^2 - 1)} (T + 1)\left[ 2T + 1 - 6T - 6h + 6\frac{(T + h)^2}{T + 1} \right] \right\} \\ & = \hat{\sigma}^2\left\{ 1 + \frac{2}{T(T - 1)} \left[ 1 - 4T - 6h + 6\frac{(T + h)^2}{(T + 1)} \right] \right\}. \end{align*}