#set document(title: "5.4 Forecast Evaluation Methods", author: "OpenStax / XYZ Homework")
#set page(width: 8.5in, height: auto, margin: 1in)
#import "@preview/cetz:0.5.2"
#set text(font: ("STIX Two Text", "Libertinus Serif", "New Computer Modern"), size: 10.5pt, lang: "en")
#show math.equation: set text(font: ("STIX Two Math", "New Computer Modern Math"))
#set par(justify: true, leading: 0.62em, spacing: 0.9em)
#set enum(spacing: 1.1em)   // room between list items so tall inline fractions don't collide
#set list(spacing: 1.1em)
#set table(stroke: 0.5pt + rgb("#c7ccd3"))

#let BLUE = rgb("#183B6F")     // brand navy — section bars + example/solution labels (white on navy 11.09:1)
#let ORANGE = rgb("#A94509")   // brand primary-700 — AA-safe deep orange for TEXT (5.93:1 on white; raw brand #F37021 is 2.94:1 and must never carry text)
#let RED = rgb("#DC2626")      // brand error-600
#let GREEN = rgb("#059669")    // brand success-600 (decoration only; small green text uses green-text #007942)

#show heading.where(level: 1): it => block(width: 100%, above: 0pt, below: 16pt,
  fill: gradient.linear(BLUE, rgb("#2C5AA0")), inset: (x: 14pt, y: 12pt), radius: 3pt,
  text(fill: white, weight: "bold", size: 19pt, it.body))
#show heading.where(level: 2): it => block(width: 100%, above: 18pt, below: 10pt,
  fill: BLUE, inset: (x: 10pt, y: 6pt), radius: 2pt,
  text(fill: white, weight: "bold", size: 12pt, it.body))
#show heading.where(level: 3): it => text(fill: ORANGE, weight: "bold", size: 12.5pt, it.body)
#show heading.where(level: 4): it => text(fill: BLUE, weight: "bold", size: 10.5pt, it.body)

#let examplebox(label, title, body) = block(width: 100%, breakable: true,
  fill: rgb("#EFF1F5"), stroke: 0.5pt + rgb("#CFDDF0"), radius: 4pt, inset: 10pt,
  above: 12pt, below: 12pt)[
  #block(below: 6pt)[#box(fill: BLUE, inset: (x: 6pt, y: 2pt), radius: 2pt,
    text(fill: white, weight: "bold", size: 8.5pt, label)) #h(0.4em) #strong[#title]]
  #body]

// rail = decorative left rule (raw brand token); labelcolor = AA-safe label text shade
#let notebox(label, rail, labelcolor, tint, body) = block(width: 100%, breakable: true, fill: tint,
    stroke: (left: 3pt + rail), inset: (left: 10pt, rest: 8pt), radius: (right: 4pt),
    above: 11pt, below: 11pt)[
    #text(fill: labelcolor, weight: "bold", size: 7.5pt, tracking: 0.5pt)[#upper(label)]
    #linebreak() #body]

#let solutionbox(body) = block(above: 4pt, below: 8pt)[
  #text(fill: BLUE, weight: "bold", size: 8.5pt)[Solution] #linebreak() #body]

#let figph(msg) = block(width: 100%, height: 60pt, fill: rgb("#f6f7f9"),
  stroke: (paint: rgb("#c7ccd3"), dash: "dashed"), radius: 4pt, inset: 10pt)[
  #align(center + horizon, text(fill: rgb("#889"), style: "italic", size: 9pt, msg))]

// Standardize inlined figure sizes: measure the natural CeTZ canvas, then scale to a
// consistent envelope (aspect-aware; see build_typst.py FIG_* constants). Unlike the
// print preamble, dimensions are FLOORED: in an editor a user can trim a figure to a
// degenerate 1-D shape (a bare line), and w/h or tw/w would then divide by zero.
#let _STD_W = 3.5
#let _WIDE_W = 5.6
#let _MAX_H = 3.4
#let _ASPECT_WIDE = 2.2
#let _UPSCALE_MAX = 1.15
#let stdfig(body) = context {
  let m = measure(body)
  let w = calc.max(m.width / 1in, 0.01)
  let h = calc.max(m.height / 1in, 0.01)
  let tw = if w / h > _ASPECT_WIDE { _WIDE_W } else { _STD_W }
  let s = calc.min(tw / w, _MAX_H / h, _UPSCALE_MAX)
  align(center, box(scale(x: s * 100%, y: s * 100%, reflow: true, body)))
}

#show figure: set block(breakable: false)
#set figure(gap: 8pt)
#show figure.caption: set text(size: 8.5pt, fill: rgb("#555"))


  
== 5.4#h(0.6em)Forecast Evaluation Methods


  
    
=== Learning Outcomes


    By the end of this section, you should be able to:


    - Explain the nature of error in forecasting a time series.
- Compute common error measures for time series models.
- Produce prediction intervals in a forecasting example.


  
  A time series forecast essentially predicts the #emph[middle of the road] values for future terms of the series. For the purposes of prediction, the model considers only the trend and any cyclical or seasonal variation that could be detected within the known data. Although noise and random variation certainly do influence future values of the time series, their precise effects cannot be predicted (by their very nature). Thus, a forecasted value is the best guess of the model in the absence of an error term. On the other hand, error does affect how certain we can be of the model’s forecasts. In this section, we focus on quantifying the error in forecasting using statistical tools such as prediction intervals.


  
    
=== Forecasting Error


    Suppose that you have a model #math.equation(block: false, alt: "open parenthesis x ^ sub n close parenthesis")[$( hat(x)_(n) )$] for a given time series #math.equation(block: false, alt: "open parenthesis x sub n close parenthesis")[$( x_(n) )$]. Recall from Components of Time Series Analysis that the residuals quantify the error between the time series the model and serve as an estimate of the noise or random variation.


    #math.equation(block: true, alt: "ε sub n equals x sub n minus x ^ sub n")[$ε_(n) = x_(n) − hat(x)_(n)$]
    The better the fit of the model to the observed data, the smaller the values of #math.equation(block: false, alt: "ε sub n")[$ε_(n)$] will be. However, even very good models can turn out to be poor predictors of future values of the time series. Ironically, the better a model does at fitting to known data, the worse it may be at predicting future data. There is a danger of #emph[overfitting]. (This term is defined in detail in Decision-Making Using Machine Learning Basics, but for now we do not need to go deeply into this topic.) A common technique to avoid overfitting is to build the model using only a portion of the known data, and then the model’s accuracy can be tested on the remaining data that was held out.


    To discuss the accuracy of a model, we should ask the complementary question: #emph[How far away are the predictions from the true values?] In other words, we should try to quantify the total error. There are several #strong[measures of error], or #strong[measures of fit], the metrics used to assess how well a model's predictions align with the observed data. We will only discuss a handful of them here that are most useful for time series. In each formula, #math.equation(block: false, alt: "x sub i")[$x_(i)$] refers to the #math.equation(block: false, alt: "i to the power th")[$i^("th")$] term of the time series, and #math.equation(block: false, alt: "ε sub i")[$ε_(i)$] is the error between the actual #math.equation(block: false, alt: "i to the power th")[$i^("th")$] term and the predicted #math.equation(block: false, alt: "i to the power th")[$i^("th")$] term of the series.


    + #strong[Mean absolute error (MAE)]: #math.equation(block: false, alt: "the fraction 1 over n ∑ i equals 1 n | ε sub i |")[$frac(1, n) limits(∑)^(n)_(i = 1) | ε_(i) |$]. A measure of the average magnitude of errors.
+ #strong[Root mean squared error (RMSE)]: #math.equation(block: false, alt: "the square root of the fraction 1 over n ∑ i equals 1 n open parenthesis ε sub i squared close parenthesis")[$sqrt(frac(1, n) limits(∑)^(n)_(i = 1) ( ε_(i)^(2) ))$]. A measure of the standard deviation of errors, penalizing larger errors more heavily than MAE does.
+ #strong[Mean absolute percentage error (MAPE)]: #math.equation(block: false, alt: "the fraction 1 over n ∑ i equals 1 n | the fraction ε sub i over x sub i |")[$frac(1, n) limits(∑)^(n)_(i = 1) | frac(ε_(i), x_(i)) |$]. A measure of the average relative errors—that is, a percentage between the predicted values and the actual values (on average).
+ #strong[Symmetric mean absolute percentage error (sMAPE)]: #math.equation(block: false, alt: "the fraction 1 over n ∑ i equals 1 n the fraction 2 | ε sub i | over | x sub i | plus | x ^ sub i |")[$frac(1, n) limits(∑)^(n)_(i = 1) frac(2 | ε_(i) |, | x_(i) | + | hat(x)_(i) |)$]. Similar to MAPE, a measure of the average relative errors, but scaled so that errors are measured in relation to both actual values and predicted values.


    Of these error measures, MAE and RMSE are #strong[scale-dependent], meaning that the error is in direct proportion to the data itself. In other words, if all terms of the time series and the model were scaled by a factor of #math.equation(block: false, alt: "k")[$k$], then the MAE and RMSE would both be multiplied by #math.equation(block: false, alt: "k")[$k$] as well. On the other hand, MAPE and sMAPE are not scale-dependent. These measures of errors are often expressed as percentages (by multiplying the result of the formula by 100%). However, neither MAPE nor sMAPE should be used for data that is measured on a scale containing 0 and negative numbers. For example, it would not be wise to use MAPE or sMAPE as a measure of error for a time series model of Celsius temperature readings. For all of these measures of error, lower values indicate less error and hence more accuracy of the model. This is useful when comparing two or more models.


    #examplebox("Example 1")[][
Compute the MAE, RMSE, MAPE, and sMAPE for the EMA smoothing model for the S&P Index time series, shown.


      
      #solutionbox[
The first step is to find the residuals, #math.equation(block: false, alt: "ε sub n")[$ε_(n)$], and take their absolute values. For RMSE, we will need the squared residuals, included in . We also include the terms #math.equation(block: false, alt: "| the fraction ε sub i over x sub i |")[$| frac(ε_(i), x_(i)) |$] and #math.equation(block: false, alt: "the fraction 2 | ε sub i | over | x sub i | plus | x ^ sub i |")[$frac(2 | ε_(i) |, | x_(i) | + | hat(x)_(i) |)$], which are used for computing MAPE and sMAPE, respectively. Notice all formulas are an average of some kind, but RMSE requires an additional step of taking a square root.


        #figure(table(
  columns: 7, align: left, inset: 6pt,
  table.header([#strong[Year]], [#strong[S&P Index at Year-End]], [#strong[EMA Estimate]], [#strong[Abs. Residuals] #math.equation(block: false, alt: "| ε sub n | equals | x sub n minus x ^ sub n |")[$attach(| ε, b: n) | = | x_(n) − hat(x)_(n) |$]], [#math.equation(block: false, alt: "open parenthesis ε sub n close parenthesis squared")[$attach(attach(( ε, b: n) ), t: 2)$]], [#math.equation(block: false, alt: "| the fraction ε sub i over x sub i |")[$| frac(ε_(i), x_(i)) |$]], [#math.equation(block: false, alt: "the fraction 2 | ε sub i | over | x sub i | plus | x ^ sub i |")[$frac(2 | ε_(i) |, | x_(i) | + | hat(x)_(i) |)$]]),
  [2013], [1848.36], [1848.36], [0], [0], [0], [0],
  [2014], [2058.9], [1848.36], [210.54], [44327.09], [0.102], [0.108],
  [2015], [2043.94], [2006.27], [37.67], [1419.03], [0.018], [0.019],
  [2016], [2238.83], [2034.52], [204.31], [41742.58], [0.091], [0.096],
  [2017], [2673.61], [2187.75], [485.86], [236059.94], [0.182], [0.2],
  [2018], [2506.85], [2552.15], [45.3], [2052.09], [0.018], [0.018],
  [2019], [3230.78], [2518.17], [712.61], [507813.01], [0.221], [0.248],
  [2020], [3756.07], [3052.63], [703.44], [494827.83], [0.187], [0.207],
  [2021], [4766.18], [3580.21], [1185.97], [1406524.84], [0.249], [0.284],
  [2022], [3839.5], [4469.69], [630.19], [397139.44], [0.164], [0.152],
  [2023], [4769.83], [3997.05], [772.78], [597188.93], [0.162], [0.176],
  [N/A], [N/A], [#strong[Average:]], [#strong[453.52]], [#strong[339008.62]], [#strong[0.127]], [#strong[0.137]],
))


        #math.equation(block: false, alt: "MAE equals 453.52")[$"MAE" #h(0.2em) = 453.52$]. (Predicted values are \[on average\] about 450 units away from true values.)


        #math.equation(block: false, alt: "RMSE equals the square root of 339008.62 equals 582.24")[$"RMSE" #h(0.2em) = sqrt(339008.62) = 582.24$]. (The errors between predicted and average values have a standard deviation of about 580 from the ideal error of 0.)


        #math.equation(block: false, alt: "MAPE equals 0.127 , or 12.7 %")[$"MAPE" #h(0.2em) = 0.127 , "or " 12.7 %$]. (Predicted values are \[on average\] within about 13% of their true values.)


        #math.equation(block: false, alt: "sMAPE equals 0.137 , or 13.7 %")[$"sMAPE" #h(0.2em) = 0.137 , "or " 13.7 %$]. (When looking at both predicted values and true values, their differences are on average about 14% from one another.)


        These results may be used to compare the EMA estimate to some other estimation/prediction method to find out which one is more accurate. The lower the errors, the more accurate the method.
]
]


  
  
    
=== Prediction Intervals


    A forecast is often accompanied by a #strong[prediction interval] giving a range of values the variable could take with some level of probability. For example, if the prediction interval of a forecast is indicated to be 80%, then the prediction interval contains a range of values that should include the actual future value with a probability of 0.8. Prediction intervals were introduced in Analysis of Variance (ANOVA) in the context of linear regression, so we won’t get into the details of the mathematics here. The key point is that we want a measure of #emph[margin of error], #math.equation(block: false, alt: "E sub n")[$E_(n)$] (depending on #math.equation(block: false, alt: "n")[$n$]), such that future observations of the data will be within the interval #math.equation(block: false, alt: "x ^ sub n plus or minus E sub n")[$hat(x)_(n) ± E_(n)$] with probability #math.equation(block: false, alt: "α")[$α$], where #math.equation(block: false, alt: "α")[$α$] is a chosen level of confidence.


    Here, we will demonstrate how to use Python to obtain prediction intervals.


    The #strong[Python] library #link("https://openstax.org/r/statsmodel")[#strong[statsmodels.tsa.arima.model]] contains functions for finding confidence intervals as well. Note that the command get\_forecast() is used here rather than forecast(), as the former contains more functionality.


    #notebox("Note", rgb("#8a94a6"), rgb("#556666"), rgb("#f7f8fa"))[
\#\#\# Please run all code from previous section before running this \#\#\#
    
    \# Set alpha to 0.2 for 80% confidence interval
    forecast\_steps = 24
    forecast\_results = results.get\_forecast(steps=forecast\_steps, alpha=0.2) 
    
    \# Extract forecast values and confidence intervals
    forecast\_values = forecast\_results.predicted\_mean
    confidence\_intervals = forecast\_results.conf\_int()
    
    \# Plot the results
    plt.figure(figsize=(10, 6))
    
    \# Plot original time series
    plt.plot(df\['Value'\], label='Original Time Series')
    
    \# Plot fitted values
    plt.plot(results.fittedvalues, color='red', label='Fitted Values')
    
    \# Plot forecasted values with confidence intervals
    plt.plot(forecast\_values, color='red', linestyle='dashed', label='Forecasted Values')
    plt.fill\_between(
      range(len(df), len(df) + forecast\_steps),
      confidence\_intervals.iloc\[:, 0\],
      confidence\_intervals.iloc\[:, 1\],
      color='red', alpha=0.2,
      label='80% Confidence Interval' 
    )
    
    \# Set labels and legend
    plt.xlabel('Months')
    plt.title('Monthly Consumption of Coal for Electricity Generation in the United States from 2016 to 2022')
    plt.legend()
    
    \# Apply the formatter to the Y-axis
    plt.gca().yaxis.set\_major\_formatter(FuncFormatter(y\_format))
    plt.show()
    
      The resulting output will look like this:


      #figure(figph[Time series plot titled Monthly consumption of coal for electricity generation in the United States from 2016 to 2022. Y-axis ranges from -50,000 to 125,000, x-axis from 0 to 100. The blue line represents the actual coal consumption, which fluctuates seasonally. The red line represents the fitted values, which smooth out the seasonal fluctuations, and the dashed red line represents the forecasted values. The shaded area represents the 80% confidence interval for the forecasted values. Coal consumption shows a general downward trend from around 125,000 tons per month in 2016 to around 0 tons per month in 2022.], alt: "Time series plot titled Monthly consumption of coal for electricity generation in the United States from 2016 to 2022. Y-axis ranges from -50,000 to 125,000, x-axis from 0 to 100. The blue line represents the actual coal consumption, which fluctuates seasonally. The red line represents the fitted values, which smooth out the seasonal fluctuations, and the dashed red line represents the forecasted values. The shaded area represents the 80% confidence interval for the forecasted values. Coal consumption shows a general downward trend from around 125,000 tons per month in 2016 to around 0 tons per month in 2022.", caption: none)
]


    The forecast data (dashed curve) is now surrounded by a shaded region. With 80% probability, all future observations should fit into the shaded region. Of course, the further into the future we try to go, the more uncertain our forecasts will be, which is indicated by the wider and wider confidence interval region.


    #notebox("Note", rgb("#8a94a6"), rgb("#556666"), rgb("#f7f8fa"))[
Datasets


      #emph[Note: The primary datasets referenced in the chapter code may also be #link("https://openstax.org/r/drive1")[downloaded here].]
]


  
  
    
      
==== Project A – Smoothing out the Stock Market


      Predicting the ups and downs of the stock market is a notoriously difficult task. On the other hand, if anyone could obtain accurate forecasts of future stock prices, then they could make a lot of money by buying up stocks at low prices and selling them later at higher prices. It is often useful just to predict whether the trend of a stock is upward or downward in the near future. In this project, you will use standard time series smoothing techniques to estimate trends in the S&P 500 Index.


      Using the data set #link("https://openstax.org/r/spreadsheetsd1j")[#strong[SP500.csv]], perform multiple techniques to smooth the time series including simple moving averages with various window lengths and exponential moving averages.


      + Produce SMA trendlines with window sizes 10, 30, 60, and 120. Graph the trendlines on the same set of axes as the original graph and describe the results. What happens to the trendline as the window size increases?
+ Produce EMA trendlines with #math.equation(block: false, alt: "α equals 0.25 , 0.5")[$α = 0.25 , 0.5$] and #math.equation(block: false, alt: "0.75")[$0.75$]. Graph the trendlines on the same set of axes as the original graph and describe the results. What happens to the trendline as the value of #math.equation(block: false, alt: "α")[$α$] increases?
+ Now compare your results from part a and part b. As a group, discuss the advantages and disadvantages of using SMA vs. EMA. What characteristics of the data are emphasized (or de-emphasized) by each method as the parameters (window size and #math.equation(block: false, alt: "α")[$α$], respectively) change.


    
    
      
==== Project B – ARIMA Analysis of Coal Consumption


      As a group, explore the time series in the data set #link("https://openstax.org/r/spreadsheetsd1j")[#strong[MonthlyCoalConsumption.xlsx]].


      + Find the smallest order of differencing so that the time series becomes stationary.
+ Then run an ACF diagram and look for any spikes that might indicate periodicity of a seasonal component.
+ Run ARIMA analysis of the time series using the parameters for #math.equation(block: false, alt: "p")[$p$] and #math.equation(block: false, alt: "d")[$d$] established by your findings about periodicity and stationarity, trying different values of #math.equation(block: false, alt: "q")[$q$]. Evaluate the model’s accuracy by computing the AIC score for each combination of parameters.
+ As a group, decide which parameters #math.equation(block: false, alt: "open parenthesis p , d , q close parenthesis")[$( p , d , q )$] you think will give the most accurate forecast. Then forecast up to 30 time units into the future and produce confidence intervals for your forecasts. Interpret your results by indicating how likely your forecasts are to be within a range of values.