(2023S) 1152-3001 Time Series Analysis

*email address: summet73y.da@gmail.com

Chapter 6. Smoothing and Decomposition methods

Section 6.1 Decomposition models

• Decomposition procedures are used in time series to describe the trend and seasonal factors in a time series.
  • More extensive decompositions might also include
    • long-run cycles,
    • holiday effects,
    • day of week effects,
    • more than one seasonal component corresponding to the different seasonal periods, and so on.
  • Here, we’ll only consider trend and seasonal decompositions.
• We will decompose a time series to three components :
  1. Trend
  2. Seasonality
  3. Random (or remainder)
• The objectives for a decomposition of a time-series are as follows :
  1. To improve understanding of the time series
  2. To improve forecast accuracy
  3. To present seasonally adjusted values by estimating seasonal effects.
     • A seasonally adjusted value removes the seasonal effect from a value so that trends can be seen more clearly.
     • For instance, in many regions of the U.S. unemployment tends to decrease in the summer due to increased employment in agricultural areas.
     • Thus a drop in the unemployment rate in June compared to May doesn’t necessarily indicate that there’s a trend toward lower unemployment in the country.
     • To see whether there is a real trend, we should adjust for the fact that unemployment is always lower in June than in May.

---

6.1.1 Basic Structure

• The following two structures are considered for basic decomposition models, where \(x_t\) is a time-series, \(S_t\) is the seasonal component, \(T_t\) is the trend component, and \(R_t\) is the random (or remainder) component
  1. Additive: \[x_t = S_t + T_t + R_t\]
  2. Multiplicative: \[x_t = S_t \times T_t \times R_t\]
• How to Choose Between Additive and Multiplicative Decompositions?
  1. The additive model : when the seasonal variation is relatively constant over time.
  2. The multiplicative model : when the variation in the seasonal pattern, or the variation around the trend, appears to be proportional to the level of the time series
     • Multiplicative decompositions are common with economic time series.
• We may consider this approach when the variations of a time series increases over time:
  1. first transform the data until the variation in the series appears to be stable over time
  2. then use an additive decomposition.
     • When a log transformation has been used, this is equivalent to using a multiplicative decomposition because \[\log{x_t} = \log{(\text{Trend})} + \log{(\text{Seasonal})} + \log{(\text{Random})}\] is equivalent to \[x_t = \text{Trend} \times \text{Seasonal} \times \text{Random}\]

---

Example 1

• In Chapter 2, we looked at quarterly beer production in Australia.
• The seasonal variation looked to be about the same magnitude across time, so an additive decomposition might be good :

beer <- scan("data//beerprod.dat") #reads the data
beer <- ts(beer)
plot(beer, type="b", ylab="beer", main="Time Series Plot of beer") # time series plot of x with points marked as "o"

---

Example 2

• The time-series of the quarterly earnings data for the Johnson and Johnson Corporations.
• The seasonal variation increases as we move across time. A multiplicative decomposition could be useful :

library(astsa)
plot(jj, type="b", ylab = 'Earnings', main = 'Quarterly Earnings of Johnson & Johnson')

---

Seasonally adjusted values

• If the seasonal component is removed from the original data, the resulting values are the “seasonally adjusted” values.
  1. Additive: the seasonally adjusted values are given by \(x_t - S_t\)
     
  2. Multiplicative: the seasonally adjusted values are given by \(x_t / S_t\)
• If the variation due to seasonality is not of primary interest, the seasonally adjusted series can be useful.
  • For instance, in many regions of the U.S. unemployment tends to decrease in the summer due to increased employment in agricultural areas.
  • Thus a drop in the unemployment rate in June compared to May doesn’t necessarily indicate that there’s a trend toward lower unemployment in the country.
  • To see whether there is a real trend, we should adjust for the fact that unemployment is always lower in June than in May.
• Seasonally adjusted series contain the random component as well as the trend component.
  • If the purpose is to interpret any changes in direction of a time-series, then it is better to use the trend-cycle component rather than the seasonally adjusted values.

---

6.1.2 Basic steps in Decomposition

1. The first step is to estimate the trend. Two different approaches could be used for this.
   • One approach is to estimate the trend with a smoothing procedure such as moving averages. With this approach, an equation is not used to describe trend.
   • The second approach is to model the trend with a regression equation.
2. The second step is to “de-trend” the series.
   • For an additive decomposition, this is done by subtracting the trend estimates from the series.
   • For a multiplicative decomposition, this is done by dividing the series by the trend values.
3. Next, seasonal components are estimated using the de-trended series.
   • For monthly data, this entails estimating an effect for each month of the year.
   • For quarterly data, this entails estimating an effect for each quarter.
   • The simplest method for estimating these effects is to average the de-trended values for a specific season.
   • For instance, to get a seasonal effect for January, we average the de-trended values for all Januarys in the series, and so on.
4. The final step is to determine the random (or remainder) component.
   • For the additive model, random = series – trend – seasonal.
   • For the multiplicative model, random = series / (trend x seasonal)

• Note that
  • the random component in the final step could be analyzed for such things as the mean location, or mean squared size (variance), or possibly even for whether the component is actually random or might be modeled with an ARIMA model.
  • A few programs iterate through the steps 1 to 3.
    • For example, after step 3 we could use the seasonal factors to de-seasonalize the series and then return to step 1 to estimate the trend based on the de-seasonalized series.

---

6.1.3 Decomposition in R

• The basic command is decompose
  • For an additive model decompose(name of series, type = "additive")
  • For a multiplicative decomposition decompose(name of series, type ="multiplicative")
• Here, the important first step is that as a preliminary, you have to use a ts command to define the seasonal span for a series.
  • For quarterly data, it might be name of series = ts(name of series, freq = 4)
  • For monthly data, it might be name of series = ts(name of series, freq = 12)
• You can plot the elements of the decomposition by putting the decompose command as an argument of a plot command.
  • As an example,

library(astsa)
decompjj <- decompose(jj, type = "multiplicative")
plot(decompjj)

• To see all elements of a stored object, simply type its name.
  • For instance, entering decompjj will show all elements of the decomposition in the example above.

decompjj

$x
          Qtr1      Qtr2      Qtr3      Qtr4
1960  0.710000  0.630000  0.850000  0.440000
1961  0.610000  0.690000  0.920000  0.550000
1962  0.720000  0.770000  0.920000  0.600000
1963  0.830000  0.800000  1.000000  0.770000
1964  0.920000  1.000000  1.240000  1.000000
1965  1.160000  1.300000  1.450000  1.250000
1966  1.260000  1.380000  1.860000  1.560000
1967  1.530000  1.590000  1.830000  1.860000
1968  1.530000  2.070000  2.340000  2.250000
1969  2.160000  2.430000  2.700000  2.250000
1970  2.790000  3.420000  3.690000  3.600000
1971  3.600000  4.320000  4.320000  4.050000
1972  4.860000  5.040000  5.040000  4.410000
1973  5.580000  5.850000  6.570000  5.310000
1974  6.030000  6.390000  6.930000  5.850000
1975  6.930000  7.740000  7.830000  6.120000
1976  7.740000  8.910000  8.280000  6.840000
1977  9.540000 10.260000  9.540000  8.729999
1978 11.880000 12.060000 12.150000  8.910000
1979 14.040000 12.960000 14.850000  9.990000
1980 16.200000 14.670000 16.020000 11.610000

$seasonal
          Qtr1      Qtr2      Qtr3      Qtr4
1960 0.9930006 1.0329845 1.1140535 0.8599614
1961 0.9930006 1.0329845 1.1140535 0.8599614
1962 0.9930006 1.0329845 1.1140535 0.8599614
1963 0.9930006 1.0329845 1.1140535 0.8599614
1964 0.9930006 1.0329845 1.1140535 0.8599614
1965 0.9930006 1.0329845 1.1140535 0.8599614
1966 0.9930006 1.0329845 1.1140535 0.8599614
1967 0.9930006 1.0329845 1.1140535 0.8599614
1968 0.9930006 1.0329845 1.1140535 0.8599614
1969 0.9930006 1.0329845 1.1140535 0.8599614
1970 0.9930006 1.0329845 1.1140535 0.8599614
1971 0.9930006 1.0329845 1.1140535 0.8599614
1972 0.9930006 1.0329845 1.1140535 0.8599614
1973 0.9930006 1.0329845 1.1140535 0.8599614
1974 0.9930006 1.0329845 1.1140535 0.8599614
1975 0.9930006 1.0329845 1.1140535 0.8599614
1976 0.9930006 1.0329845 1.1140535 0.8599614
1977 0.9930006 1.0329845 1.1140535 0.8599614
1978 0.9930006 1.0329845 1.1140535 0.8599614
1979 0.9930006 1.0329845 1.1140535 0.8599614
1980 0.9930006 1.0329845 1.1140535 0.8599614

$trend
         Qtr1     Qtr2     Qtr3     Qtr4
1960       NA       NA  0.64500  0.64000
1961  0.65625  0.67875  0.70625  0.73000
1962  0.74000  0.74625  0.76625  0.78375
1963  0.79750  0.82875  0.86125  0.89750
1964  0.95250  1.01125  1.07000  1.13750
1965  1.20125  1.25875  1.30250  1.32500
1966  1.38625  1.47625  1.54875  1.60875
1967  1.63125  1.66500  1.70250  1.76250
1968  1.88625  1.99875  2.12625  2.25000
1969  2.34000  2.38500  2.46375  2.66625
1970  2.91375  3.20625  3.47625  3.69000
1971  3.88125  4.01625  4.23000  4.47750
1972  4.65750  4.79250  4.92750  5.11875
1973  5.41125  5.71500  5.88375  6.00750
1974  6.12000  6.23250  6.41250  6.69375
1975  6.97500  7.12125  7.25625  7.50375
1976  7.70625  7.85250  8.16750  8.56125
1977  8.88750  9.28125  9.81000 10.32750
1978 10.87875 11.22750 11.52000 11.90250
1979 12.35250 12.82500 13.23000 13.71375
1980 14.07375 14.42250       NA       NA

$random
          Qtr1      Qtr2      Qtr3      Qtr4
1960        NA        NA 1.1829139 0.7994545
1961 0.9360758 0.9841141 1.1692929 0.8761145
1962 0.9798312 0.9988783 1.0777332 0.8902147
1963 1.0480883 0.9344857 1.0422327 0.9976480
1964 0.9726875 0.9572991 1.0402359 1.0222795
1965 0.9724675 0.9997929 0.9992731 1.0970216
1966 0.9153338 0.9049516 1.0780169 1.1276053
1967 0.9445423 0.9244620 0.9648458 1.2271704
1968 0.8168507 1.0025778 0.9878602 1.1628429
1969 0.9295835 0.9863342 0.9836964 0.9813020
1970 0.9642783 1.0326067 0.9528166 1.1344809
1971 0.9340742 1.0412840 0.9167213 1.0518177
1972 1.0508335 1.0180629 0.9181166 1.0018339
1973 1.0384536 0.9909365 1.0023170 1.0278312
1974 0.9922392 0.9925326 0.9700627 1.0162661
1975 1.0005517 1.0521821 0.9685978 0.9484056
1976 1.0114592 1.0984390 0.9099869 0.9290519
1977 1.0809840 1.0701560 0.8729177 0.9829695
1978 1.0997347 1.0398494 0.9467117 0.8704835
1979 1.1446237 0.9782589 1.0075359 0.8470915
1980 1.1591928 0.9846815        NA        NA

$figure
[1] 0.9930006 1.0329845 1.1140535 0.8599614

$type
[1] "multiplicative"

attr(,"class")
[1] "decomposed.ts"

• When the decomposition is stored in an object, you also have access to the various elements of the decomposition.
  • For instance, decompjj$figure contains the seasonal effects values for four quarters :

decompjj$figure

[1] 0.9930006 1.0329845 1.1140535 0.8599614

---

Example 1 (Continued)

• The following commands produced the graph and numerical output that follows for the Australian beer production series.

beer <- scan("data//beerprod.dat")
beer <- ts(beer, freq = 4)
decompbeer <- decompose(beer, type = "additive")
plot(decompbeer)

decompbeer

$x
    Qtr1  Qtr2  Qtr3  Qtr4
1  284.4 212.8 226.9 308.4
2  262.0 227.9 236.1 320.4
3  271.9 232.8 237.0 313.4
4  261.4 226.8 249.9 314.3
5  286.1 226.5 260.4 311.4
6  294.7 232.6 257.2 339.2
7  279.1 249.8 269.8 345.7
8  293.8 254.7 277.5 363.4
9  313.4 272.8 300.1 369.5
10 330.8 287.8 305.9 386.1
11 335.2 288.0 308.3 402.3
12 352.8 316.1 324.9 404.8
13 393.0 318.9 327.0 442.3
14 383.1 331.6 361.4 445.9
15 386.6 357.2 373.6 466.2
16 409.6 369.8 378.6 487.0
17 419.2 376.7 392.8 506.1
18 458.4 387.4 426.9 525.0

$seasonal
         Qtr1       Qtr2       Qtr3       Qtr4
1    7.896324 -40.678676 -24.650735  57.433088
2    7.896324 -40.678676 -24.650735  57.433088
3    7.896324 -40.678676 -24.650735  57.433088
4    7.896324 -40.678676 -24.650735  57.433088
5    7.896324 -40.678676 -24.650735  57.433088
6    7.896324 -40.678676 -24.650735  57.433088
7    7.896324 -40.678676 -24.650735  57.433088
8    7.896324 -40.678676 -24.650735  57.433088
9    7.896324 -40.678676 -24.650735  57.433088
10   7.896324 -40.678676 -24.650735  57.433088
11   7.896324 -40.678676 -24.650735  57.433088
12   7.896324 -40.678676 -24.650735  57.433088
13   7.896324 -40.678676 -24.650735  57.433088
14   7.896324 -40.678676 -24.650735  57.433088
15   7.896324 -40.678676 -24.650735  57.433088
16   7.896324 -40.678676 -24.650735  57.433088
17   7.896324 -40.678676 -24.650735  57.433088
18   7.896324 -40.678676 -24.650735  57.433088

$trend
       Qtr1     Qtr2     Qtr3     Qtr4
1        NA       NA 255.3250 254.4125
2  257.4500 260.1000 262.8375 264.6875
3  265.4125 264.6500 262.4625 260.4000
4  261.2625 262.9875 266.1875 269.2375
5  270.5125 271.4625 272.1750 274.0125
6  274.3750 277.4500 278.9750 279.1750
7  282.9000 285.2875 287.9375 290.3875
8  291.9625 295.1375 299.8000 304.5125
9  309.6000 313.1875 316.1250 320.1750
10 322.7750 325.5750 328.2000 328.7750
11 329.1000 331.4250 335.6500 341.3625
12 346.9500 349.3375 354.6750 360.0500
13 360.6625 365.6125 369.0625 369.4125
14 375.3000 380.0500 380.9375 384.5750
15 389.3000 393.3625 398.7750 403.2250
16 405.4250 408.6500 412.4500 414.5125
17 417.1500 421.3125 428.6000 434.8375
18 440.4375 447.0625       NA       NA

$random
           Qtr1         Qtr2         Qtr3         Qtr4
1            NA           NA  -3.77426471  -3.44558824
2   -3.34632353   8.47867647  -2.08676471  -1.72058824
3   -1.40882353   8.82867647  -0.81176471  -4.43308824
4   -7.75882353   4.49117647   8.36323529 -12.37058824
5    7.69117647  -4.28382353  12.87573529 -20.04558824
6   12.42867647  -4.17132353   2.87573529   2.59191176
7  -11.69632353   5.19117647   6.51323529  -2.12058824
8   -6.05882353   0.24117647   2.35073529   1.45441176
9   -4.09632353   0.29117647   8.62573529  -8.10808824
10   0.12867647   2.90367647   2.35073529  -0.10808824
11  -1.79632353  -2.74632353  -2.69926471   3.50441176
12  -2.04632353   7.44117647  -5.12426471 -12.68308824
13  24.44117647  -6.03382353 -17.41176471  15.45441176
14  -0.09632353  -7.77132353   5.11323529   3.89191176
15 -10.59632353   4.51617647  -0.52426471   5.54191176
16  -3.72132353   1.82867647  -9.19926471  15.05441176
17  -5.84632353  -3.93382353 -11.14926471  13.82941176
18  10.06617647 -18.98382353           NA           NA

$figure
[1]   7.896324 -40.678676 -24.650735  57.433088

$type
[1] "additive"

attr(,"class")
[1] "decomposed.ts"

• The plot shows

  1. the observed series,
  2. the smoothed trend line,
  3. the seasonal pattern,
  4. the random part of the series.

• Note that the elements of $figure are the effects for the four quarters. Then the seasonal effect values are repeated each year (row) in the $seasonal object at the top of this page.

• The seasonal values can be used to seasonally adjusted future values.

  • Suppose for example that the next quarter 4 seasonal value past the end of the series has the value 535.
  • The quarter 4 seasonal effect is about 57.43.
  • Thus for this future value, the seasonally adjusted value will be \(535 − 57.43 = 477.57\).

---

Example 2 (Continued)

• The following two commands will do
  • a multiplicative decomposition of the quarterly earnings data for the Johnson and Johnson Corporations and
  • print the seasonal effects.

decombeermult <- decompose(jj, type = "multiplicative")
decombeermult$figure

[1] 0.9930006 1.0329845 1.1140535 0.8599614

• To seasonally adjust a value, divide the observed value of the series by the seasonal factors.
  • For example if a future quarter 4 value is 10, the seasonally adjusted value = 10/0.8599614 = 11.62843.

---

How the Trend Values Were Calculated : Moving average

• The classical method of time series decomposition uses a moving average method to estimate the trend component :
• Moving average of order \(m=2k+1\) : \[ \hat{T}_t = {1 \over m} \sum_{=-k}^{k}{x_{t+j}} \]
  • (Example) For a daily time-series data with a weekly seasonality (Frequency = 7), moving average of order \(m=7\) can be used : \[ \hat{T}_t = {1 \over 7}(x_{t-3} + x_{t-2} + \dots + x_{t+3}) \]
• Moving average of order \(m=2k\) : We apply a moving average to a moving average \[ \hat{T}_t = {1 \over 2m} x_{t-k} + {1 \over m} \sum_{-k+1}^{k-1}{x_{t+j}} + {1 \over 2m} x_{t+k} \]
  • (Example) For a quarterly time-series data (Frequency = 4), we have the following :
    • Here’s how the centered moving average for \(t = 3\) would be calculated.
    • Average the observed data values at times 1 to 4: \[ {1\over 4}(x_1 + x_2 + x_3 + x_4) \]
    • Then average those two averages: \[\begin{equation} \begin{split} & {1\over 2} \left( {1\over 4}(x_1 + x_2 + x_3 + x_4) + {1\over 4}(x_2 + x_3 + x_4 + x_5) \right) \\ & = {1 \over 8}x_1 + {1 \over 4}x_2 + {1 \over 4}x_3 + {1 \over 4}x_4 + {1 \over 8}x_5 \end{split} \end{equation}\]
  • More generally, the moving average for time \(t\) (with 4 quarters) is \[ {1 \over 8}x_{t-2} + {1 \over 4}x_{t-1} + {1 \over 4}x_t + {1 \over 4}x_{t+1} + {1 \over 8}x_{t+2} \]
  • For monthly data the moving average for time \(t\) will be \[ {1 \over 24}x_{t-6} + {1 \over 12}\sum_{j=t-5}^{t+5}{x_{j}} + {1 \over 24}x_{t+6} \]

---

Decomposition with moving average

• Additive decomposition
  1. Estimate the trend : compute the trend component \(\hat{T}_t\)
     • If the frequency of the time-series is odd, obtain the moving average of order \(m=2k+1\) : \[ \hat{T}_t = {1 \over m} \sum_{=-k}^{k}{x_{t+j}} \]
     • If the frequency of the time-series is even, obtain the moving average of order \(m=2k\) : \[ \hat{T}_t = {1 \over 2m} x_{t-k} + {1 \over m} \sum_{-k+1}^{k-1}{x_{t+j}} + {1 \over 2m} x_{t+k} \]
  2. Calculate the de-trended series : \(x_t - \hat{T}_t\)
  3. Estimate the seasonal components for each season : \(\hat{S}_t\)
     • Simply average the de-trended values for each season.
     • These seasonal component values are then adjusted by adding a constant to ensure that they add to zero
  4. Determine the random component : \(\hat{R}_t = x_t - \hat{T}_t - \hat{S}_t\)

### library
library(astsa)
library(dplyr)

### Example 1
beer <- scan("data//beerprod.dat")
beer <- ts(beer, freq = 4)
decompbeer <- decompose(beer, type = "additive")
plot(decompbeer)

# trend
beer %>% stats::filter(filter = c(1/8, 1/4, 1/4, 1/4, 1/8), sides=2) %>% as.numeric()

 [1]       NA       NA 255.3250 254.4125 257.4500 260.1000 262.8375
 [8] 264.6875 265.4125 264.6500 262.4625 260.4000 261.2625 262.9875
[15] 266.1875 269.2375 270.5125 271.4625 272.1750 274.0125 274.3750
[22] 277.4500 278.9750 279.1750 282.9000 285.2875 287.9375 290.3875
[29] 291.9625 295.1375 299.8000 304.5125 309.6000 313.1875 316.1250
[36] 320.1750 322.7750 325.5750 328.2000 328.7750 329.1000 331.4250
[43] 335.6500 341.3625 346.9500 349.3375 354.6750 360.0500 360.6625
[50] 365.6125 369.0625 369.4125 375.3000 380.0500 380.9375 384.5750
[57] 389.3000 393.3625 398.7750 403.2250 405.4250 408.6500 412.4500
[64] 414.5125 417.1500 421.3125 428.6000 434.8375 440.4375 447.0625
[71]       NA       NA

decompbeer$trend %>% as.numeric()

 [1]       NA       NA 255.3250 254.4125 257.4500 260.1000 262.8375
 [8] 264.6875 265.4125 264.6500 262.4625 260.4000 261.2625 262.9875
[15] 266.1875 269.2375 270.5125 271.4625 272.1750 274.0125 274.3750
[22] 277.4500 278.9750 279.1750 282.9000 285.2875 287.9375 290.3875
[29] 291.9625 295.1375 299.8000 304.5125 309.6000 313.1875 316.1250
[36] 320.1750 322.7750 325.5750 328.2000 328.7750 329.1000 331.4250
[43] 335.6500 341.3625 346.9500 349.3375 354.6750 360.0500 360.6625
[50] 365.6125 369.0625 369.4125 375.3000 380.0500 380.9375 384.5750
[57] 389.3000 393.3625 398.7750 403.2250 405.4250 408.6500 412.4500
[64] 414.5125 417.1500 421.3125 428.6000 434.8375 440.4375 447.0625
[71]       NA       NA

# seasonality
decompbeer$seasonal[1:4]

[1]   7.896324 -40.678676 -24.650735  57.433088

decompbeer$seasonal[1:4] %>% sum()

[1] 3.552714e-15

tmp.matrix <- matrix( decompbeer$x - decompbeer$trend, ncol=4, byrow=T )
tmp.matrix %>% colMeans(na.rm=T)

[1]   7.677941 -40.897059 -24.869118  57.214706

tmp.matrix %>% colMeans(na.rm=T) %>% sum()

[1] -0.8735294

decompbeer$seasonal[1:4] - (tmp.matrix %>% colMeans(na.rm=T))

[1] 0.2183824 0.2183824 0.2183824 0.2183824

# remainder
decompbeer$x - decompbeer$trend - decompbeer$seasonal

           Qtr1         Qtr2         Qtr3         Qtr4
1            NA           NA  -3.77426471  -3.44558824
2   -3.34632353   8.47867647  -2.08676471  -1.72058824
3   -1.40882353   8.82867647  -0.81176471  -4.43308824
4   -7.75882353   4.49117647   8.36323529 -12.37058824
5    7.69117647  -4.28382353  12.87573529 -20.04558824
6   12.42867647  -4.17132353   2.87573529   2.59191176
7  -11.69632353   5.19117647   6.51323529  -2.12058824
8   -6.05882353   0.24117647   2.35073529   1.45441176
9   -4.09632353   0.29117647   8.62573529  -8.10808824
10   0.12867647   2.90367647   2.35073529  -0.10808824
11  -1.79632353  -2.74632353  -2.69926471   3.50441176
12  -2.04632353   7.44117647  -5.12426471 -12.68308824
13  24.44117647  -6.03382353 -17.41176471  15.45441176
14  -0.09632353  -7.77132353   5.11323529   3.89191176
15 -10.59632353   4.51617647  -0.52426471   5.54191176
16  -3.72132353   1.82867647  -9.19926471  15.05441176
17  -5.84632353  -3.93382353 -11.14926471  13.82941176
18  10.06617647 -18.98382353           NA           NA

decompbeer$random

           Qtr1         Qtr2         Qtr3         Qtr4
1            NA           NA  -3.77426471  -3.44558824
2   -3.34632353   8.47867647  -2.08676471  -1.72058824
3   -1.40882353   8.82867647  -0.81176471  -4.43308824
4   -7.75882353   4.49117647   8.36323529 -12.37058824
5    7.69117647  -4.28382353  12.87573529 -20.04558824
6   12.42867647  -4.17132353   2.87573529   2.59191176
7  -11.69632353   5.19117647   6.51323529  -2.12058824
8   -6.05882353   0.24117647   2.35073529   1.45441176
9   -4.09632353   0.29117647   8.62573529  -8.10808824
10   0.12867647   2.90367647   2.35073529  -0.10808824
11  -1.79632353  -2.74632353  -2.69926471   3.50441176
12  -2.04632353   7.44117647  -5.12426471 -12.68308824
13  24.44117647  -6.03382353 -17.41176471  15.45441176
14  -0.09632353  -7.77132353   5.11323529   3.89191176
15 -10.59632353   4.51617647  -0.52426471   5.54191176
16  -3.72132353   1.82867647  -9.19926471  15.05441176
17  -5.84632353  -3.93382353 -11.14926471  13.82941176
18  10.06617647 -18.98382353           NA           NA

• Multiplicative decomposition
  1. Estimate the trend : compute the trend component \(\hat{T}_t\)
     • If the frequency of the time-series is odd, obtain the moving average of order \(m=2k+1\) : \[ \hat{T}_t = {1 \over m} \sum_{=-k}^{k}{x_{t+j}} \]
     • If the frequency of the time-series is even, obtain the moving average of order \(m=2k\) : \[ \hat{T}_t = {1 \over 2m} x_{t-k} + {1 \over m} \sum_{-k+1}^{k-1}{x_{t+j}} + {1 \over 2m} x_{t+k} \]
  2. Calculate the de-trended series : \(x_t / \hat{T}_t\)
  3. Estimate the seasonal components for each season : \(\hat{S}_t\)
     • Simply average the de-trended values for each season.
     • These seasonal component values are then adjusted by multiplying a constant to ensure that they add to \(m\)
  4. Determine the random component : \(\hat{R}_t = x_t / (\hat{T}_t \times \hat{S}_t)\)

### Example 2
decompjj <- decompose(jj, type = "multiplicative")
plot(decompjj)

# trend
jj %>% stats::filter(filter = c(1/8, 1/4, 1/4, 1/4, 1/8), sides=2) %>% as.numeric()

 [1]       NA       NA  0.64500  0.64000  0.65625  0.67875  0.70625
 [8]  0.73000  0.74000  0.74625  0.76625  0.78375  0.79750  0.82875
[15]  0.86125  0.89750  0.95250  1.01125  1.07000  1.13750  1.20125
[22]  1.25875  1.30250  1.32500  1.38625  1.47625  1.54875  1.60875
[29]  1.63125  1.66500  1.70250  1.76250  1.88625  1.99875  2.12625
[36]  2.25000  2.34000  2.38500  2.46375  2.66625  2.91375  3.20625
[43]  3.47625  3.69000  3.88125  4.01625  4.23000  4.47750  4.65750
[50]  4.79250  4.92750  5.11875  5.41125  5.71500  5.88375  6.00750
[57]  6.12000  6.23250  6.41250  6.69375  6.97500  7.12125  7.25625
[64]  7.50375  7.70625  7.85250  8.16750  8.56125  8.88750  9.28125
[71]  9.81000 10.32750 10.87875 11.22750 11.52000 11.90250 12.35250
[78] 12.82500 13.23000 13.71375 14.07375 14.42250       NA       NA

decompjj$trend %>% as.numeric()

 [1]       NA       NA  0.64500  0.64000  0.65625  0.67875  0.70625
 [8]  0.73000  0.74000  0.74625  0.76625  0.78375  0.79750  0.82875
[15]  0.86125  0.89750  0.95250  1.01125  1.07000  1.13750  1.20125
[22]  1.25875  1.30250  1.32500  1.38625  1.47625  1.54875  1.60875
[29]  1.63125  1.66500  1.70250  1.76250  1.88625  1.99875  2.12625
[36]  2.25000  2.34000  2.38500  2.46375  2.66625  2.91375  3.20625
[43]  3.47625  3.69000  3.88125  4.01625  4.23000  4.47750  4.65750
[50]  4.79250  4.92750  5.11875  5.41125  5.71500  5.88375  6.00750
[57]  6.12000  6.23250  6.41250  6.69375  6.97500  7.12125  7.25625
[64]  7.50375  7.70625  7.85250  8.16750  8.56125  8.88750  9.28125
[71]  9.81000 10.32750 10.87875 11.22750 11.52000 11.90250 12.35250
[78] 12.82500 13.23000 13.71375 14.07375 14.42250       NA       NA

# seasonality
decompjj$seasonal[1:4]

[1] 0.9930006 1.0329845 1.1140535 0.8599614

decompjj$seasonal[1:4] %>% sum()

[1] 4

tmp.matrix <- matrix( decompjj$x / decompjj$trend, ncol=4, byrow=T )
tmp.matrix %>% colMeans(na.rm=T)

[1] 0.9925977 1.0325654 1.1136015 0.8596125

tmp.matrix %>% colMeans(na.rm=T) %>% sum()

[1] 3.998377

decompjj$seasonal[1:4] / (tmp.matrix %>% colMeans(na.rm=T))

[1] 1.000406 1.000406 1.000406 1.000406

# remainder
decompjj$x / decompjj$trend / decompjj$seasonal

          Qtr1      Qtr2      Qtr3      Qtr4
1960        NA        NA 1.1829139 0.7994545
1961 0.9360758 0.9841141 1.1692929 0.8761145
1962 0.9798312 0.9988783 1.0777332 0.8902147
1963 1.0480883 0.9344857 1.0422327 0.9976480
1964 0.9726875 0.9572991 1.0402359 1.0222795
1965 0.9724675 0.9997929 0.9992731 1.0970216
1966 0.9153338 0.9049516 1.0780169 1.1276053
1967 0.9445423 0.9244620 0.9648458 1.2271704
1968 0.8168507 1.0025778 0.9878602 1.1628429
1969 0.9295835 0.9863342 0.9836964 0.9813020
1970 0.9642783 1.0326067 0.9528166 1.1344809
1971 0.9340742 1.0412840 0.9167213 1.0518177
1972 1.0508335 1.0180629 0.9181166 1.0018339
1973 1.0384536 0.9909365 1.0023170 1.0278312
1974 0.9922392 0.9925326 0.9700627 1.0162661
1975 1.0005517 1.0521821 0.9685978 0.9484056
1976 1.0114592 1.0984390 0.9099869 0.9290519
1977 1.0809840 1.0701560 0.8729177 0.9829695
1978 1.0997347 1.0398494 0.9467117 0.8704835
1979 1.1446237 0.9782589 1.0075359 0.8470915
1980 1.1591928 0.9846815        NA        NA

decompjj$random

          Qtr1      Qtr2      Qtr3      Qtr4
1960        NA        NA 1.1829139 0.7994545
1961 0.9360758 0.9841141 1.1692929 0.8761145
1962 0.9798312 0.9988783 1.0777332 0.8902147
1963 1.0480883 0.9344857 1.0422327 0.9976480
1964 0.9726875 0.9572991 1.0402359 1.0222795
1965 0.9724675 0.9997929 0.9992731 1.0970216
1966 0.9153338 0.9049516 1.0780169 1.1276053
1967 0.9445423 0.9244620 0.9648458 1.2271704
1968 0.8168507 1.0025778 0.9878602 1.1628429
1969 0.9295835 0.9863342 0.9836964 0.9813020
1970 0.9642783 1.0326067 0.9528166 1.1344809
1971 0.9340742 1.0412840 0.9167213 1.0518177
1972 1.0508335 1.0180629 0.9181166 1.0018339
1973 1.0384536 0.9909365 1.0023170 1.0278312
1974 0.9922392 0.9925326 0.9700627 1.0162661
1975 1.0005517 1.0521821 0.9685978 0.9484056
1976 1.0114592 1.0984390 0.9099869 0.9290519
1977 1.0809840 1.0701560 0.8729177 0.9829695
1978 1.0997347 1.0398494 0.9467117 0.8704835
1979 1.1446237 0.9782589 1.0075359 0.8470915
1980 1.1591928 0.9846815        NA        NA

### Example 2-2
decompjj2 <- decompose(log(jj), type = "additive")

# trend
log(jj) %>% stats::filter(filter = c(1/8, 1/4, 1/4, 1/4, 1/8), sides=2) %>% as.numeric()

 [1]            NA            NA -0.4659820641 -0.4735863435
 [5] -0.4523227061 -0.4145375971 -0.3659206214 -0.3314842249
 [9] -0.3177718602 -0.3068954381 -0.2782472049 -0.2556977422
[13] -0.2404973894 -0.1988920809 -0.1548409773 -0.1140795372
[17] -0.0592976708  0.0002618472  0.0619076444  0.1236783792
[21]  0.1760304344  0.2234799004  0.2617093088  0.2795106776
[25]  0.3181021984  0.3769215985  0.4288838841  0.4708597005
[29]  0.4865334500  0.5064872182  0.5284735516  0.5614503754
[33]  0.6251565696  0.6796801560  0.7465794329  0.8097273249
[37]  0.8476577616  0.8655453671  0.8975370388  0.9722473723
[41]  1.0540128697  1.1518101589  1.2424221438  1.3034855314
[45]  1.3523910058  1.3868175033  1.4390534568  1.4958353659
[49]  1.5343730358  1.5642865968  1.5922001152  1.6280983549
[53]  1.6798662711  1.7362193831  1.7691288056  1.7898601608
[57]  1.8075649843  1.8263394601  1.8558347887  1.8971823164
[61]  1.9364035809  1.9573064724  1.9767642611  2.0081785647
[65]  2.0327601912  2.0536484529  2.0936881320  2.1374594315
[69]  2.1728005709  2.2210040760  2.2789216201  2.3265471426
[73]  2.3769817720  2.4097623559  2.4331952397  2.4630736877
[77]  2.4971542121  2.5365393430  2.5687282423  2.6021079604
[81]  2.6270798326  2.6553448676            NA            NA

decompjj2$trend %>% as.numeric()

 [1]            NA            NA -0.4659820641 -0.4735863435
 [5] -0.4523227061 -0.4145375971 -0.3659206214 -0.3314842249
 [9] -0.3177718602 -0.3068954381 -0.2782472049 -0.2556977422
[13] -0.2404973894 -0.1988920809 -0.1548409773 -0.1140795372
[17] -0.0592976708  0.0002618472  0.0619076444  0.1236783792
[21]  0.1760304344  0.2234799004  0.2617093088  0.2795106776
[25]  0.3181021984  0.3769215985  0.4288838841  0.4708597005
[29]  0.4865334500  0.5064872182  0.5284735516  0.5614503754
[33]  0.6251565696  0.6796801560  0.7465794329  0.8097273249
[37]  0.8476577616  0.8655453671  0.8975370388  0.9722473723
[41]  1.0540128697  1.1518101589  1.2424221438  1.3034855314
[45]  1.3523910058  1.3868175033  1.4390534568  1.4958353659
[49]  1.5343730358  1.5642865968  1.5922001152  1.6280983549
[53]  1.6798662711  1.7362193831  1.7691288056  1.7898601608
[57]  1.8075649843  1.8263394601  1.8558347887  1.8971823164
[61]  1.9364035809  1.9573064724  1.9767642611  2.0081785647
[65]  2.0327601912  2.0536484529  2.0936881320  2.1374594315
[69]  2.1728005709  2.2210040760  2.2789216201  2.3265471426
[73]  2.3769817720  2.4097623559  2.4331952397  2.4630736877
[77]  2.4971542121  2.5365393430  2.5687282423  2.6021079604
[81]  2.6270798326  2.6553448676            NA            NA

# check
decompjj2$trend %>% as.numeric() %>% exp()

 [1]         NA         NA  0.6275185  0.6227648  0.6361488  0.6606457
 [7]  0.6935579  0.7178575  0.7277688  0.7357275  0.7571096  0.7743760
[13]  0.7862367  0.8196383  0.8565514  0.8921870  0.9424262  1.0002619
[19]  1.0638641  1.1316518  1.1924744  1.2504205  1.2991488  1.3224825
[25]  1.3745167  1.4577900  1.5355427  1.6013703  1.6266675  1.6594517
[31]  1.6963410  1.7532135  1.8685385  1.9732465  2.1097710  2.2472951
[37]  2.3341733  2.3763017  2.4535527  2.6438796  2.8691415  3.1639149
[43]  3.4639936  3.6821084  3.8666597  4.0020931  4.2167026  4.4630633
[49]  4.6384165  4.7792642  4.9145496  5.0941782  5.3648385  5.6758446
[55]  5.8657409  5.9886150  6.0955865  6.2111090  6.3970362  6.6670822
[61]  6.9337693  7.0802306  7.2193452  7.4497358  7.6351317  7.7962937
[67]  8.1147885  8.4778716  8.7828466  9.2165804  9.7661431 10.2425145
[73] 10.7723404 11.1313155 11.3952345 11.7408438 12.1478745 12.6358668
[79] 13.0492184 13.4921490 13.8333153 14.2298926         NA         NA

decompjj$trend %>% as.numeric()

 [1]       NA       NA  0.64500  0.64000  0.65625  0.67875  0.70625
 [8]  0.73000  0.74000  0.74625  0.76625  0.78375  0.79750  0.82875
[15]  0.86125  0.89750  0.95250  1.01125  1.07000  1.13750  1.20125
[22]  1.25875  1.30250  1.32500  1.38625  1.47625  1.54875  1.60875
[29]  1.63125  1.66500  1.70250  1.76250  1.88625  1.99875  2.12625
[36]  2.25000  2.34000  2.38500  2.46375  2.66625  2.91375  3.20625
[43]  3.47625  3.69000  3.88125  4.01625  4.23000  4.47750  4.65750
[50]  4.79250  4.92750  5.11875  5.41125  5.71500  5.88375  6.00750
[57]  6.12000  6.23250  6.41250  6.69375  6.97500  7.12125  7.25625
[64]  7.50375  7.70625  7.85250  8.16750  8.56125  8.88750  9.28125
[71]  9.81000 10.32750 10.87875 11.22750 11.52000 11.90250 12.35250
[78] 12.82500 13.23000 13.71375 14.07375 14.42250       NA       NA

decompjj2$seasonal[1:4] %>% as.numeric() %>% exp()

[1] 0.9990446 1.0404210 1.1178227 0.8606629

decompjj$seasonal[1:4] %>% as.numeric()

[1] 0.9930006 1.0329845 1.1140535 0.8599614

decompjj2$random %>% as.numeric() %>% exp()

 [1]        NA        NA 1.2117678 0.8209100 0.9598121 1.0038559
 [7] 1.1866761 0.8902078 0.9902713 1.0059227 1.0870665 0.9002566
[13] 1.0566713 0.9381205 1.0444163 1.0027710 0.9771373 0.9608977
[19] 1.0427077 1.0267249 0.9736975 0.9992592 0.9984726 1.0982140
[25] 0.9175625 0.9098609 1.0836227 1.1318784 0.9414728 0.9209232
[31] 0.9650837 1.2326650 0.8196048 1.0082771 0.9922191 1.1632936
[37] 0.9262661 0.9828689 0.9844540 0.9887984 0.9733462 1.0389441
[43] 0.9529637 1.1359857 0.9319265 1.0374984 0.9165113 1.0543600
[49] 1.0487734 1.0135855 0.9174320 1.0058457 1.0411005 0.9906409
[55] 1.0020042 1.0302320 0.9901864 0.9888322 0.9691288 1.0194994
[61] 1.0004122 1.0507138 0.9702664 0.9545034 1.0147044 1.0984503
[67] 0.9128096 0.9374243 1.0872470 1.0699622 0.8738811 0.9903177
[73] 1.1038791 1.0413380 0.9538500 0.8817497 1.1568630 0.9858046
[79] 1.0180498 0.8603027 1.1722058 0.9908762        NA        NA

decompjj$random %>% as.numeric()

 [1]        NA        NA 1.1829139 0.7994545 0.9360758 0.9841141
 [7] 1.1692929 0.8761145 0.9798312 0.9988783 1.0777332 0.8902147
[13] 1.0480883 0.9344857 1.0422327 0.9976480 0.9726875 0.9572991
[19] 1.0402359 1.0222795 0.9724675 0.9997929 0.9992731 1.0970216
[25] 0.9153338 0.9049516 1.0780169 1.1276053 0.9445423 0.9244620
[31] 0.9648458 1.2271704 0.8168507 1.0025778 0.9878602 1.1628429
[37] 0.9295835 0.9863342 0.9836964 0.9813020 0.9642783 1.0326067
[43] 0.9528166 1.1344809 0.9340742 1.0412840 0.9167213 1.0518177
[49] 1.0508335 1.0180629 0.9181166 1.0018339 1.0384536 0.9909365
[55] 1.0023170 1.0278312 0.9922392 0.9925326 0.9700627 1.0162661
[61] 1.0005517 1.0521821 0.9685978 0.9484056 1.0114592 1.0984390
[67] 0.9099869 0.9290519 1.0809840 1.0701560 0.8729177 0.9829695
[73] 1.0997347 1.0398494 0.9467117 0.8704835 1.1446237 0.9782589
[79] 1.0075359 0.8470915 1.1591928 0.9846815        NA        NA

---

Comments on the decomposition method with moving average

• While the classical decomposition (using moving average for trend) is still widely used, it is not recommended, as there are now several much better methods.
• Some of the problems with classical decomposition are as follows:
  • The estimate of the trend component is unavailable for the first few and last few observations.
  • The trend components estimate tends to over-smooth rapid rises and falls in the data
  • The classical decomposition method assume that the seasonal component repeats from year to year.
    • For many series, this is a reasonable assumption, but for some longer series it is not.
    • For example, electricity demand patterns have changed over time as air conditioning has become more widespread.
      • Specifically, in many locations, the seasonal usage pattern from several decades ago had its maximum demand in winter (due to heating), while the current seasonal pattern has its maximum demand in summer (due to air conditioning).
    • The classical decomposition methods are unable to capture these seasonal changes over time.
  • The classical decomposition method is not robust to outliers
    • For example, the monthly air passenger traffic may be affected by an industrial dispute, making the traffic during the dispute different from usual.

---

6.1.5 STL decomposition

• STL is a versatile and robust method for decomposing time series.
• STL is an acronym for “Seasonal and Trend decomposition using Loess”
  • Loess means the locally estimated scatterplot smoothing.
  • Loess is one approach for local regression
• STL has several advantages :
  • The seasonal component is allowed to change over time, and the rate of change can be controlled by the user.
  • The smoothness of the trend component can also be controlled by the user.
  • It can be robust to outliers (i.e., the user can specify a robust decomposition)
    • Outliers will not affect the estimates of the trend and seasonal components. They will, however, affect the remainder component.
• The R command stl() does an additive decomposition in which a loess smoother is used to estimate the trend and (potentially) the seasonal effects as well.
  • There are several parameters that can be adjusted, but the default does a fairly good job.
  • It is possible to obtain a multiplicative decomposition by first taking logs of the data, then back-transforming the components.
  • Decompositions between additive and multiplicative can be obtained using a Box-Cox transformation of the data.

Box-Cox transformation

• If a time-series shows variation that increases or decreases with the level of the series, then a transformation can be useful.
  • (Example) logarithmic transformation : consider \(w_t = \log(x_t)\)
    • Logarithms are useful because they are interpretable : changes in a log value are relative (or percentage) changes on the original scale.
  • (Example) power transformations : consider \(w_t = x_t^p\) with a constant \(p\) (like \(p=1/2\) or \(2\))
• A useful family of transformations, that includes both logarithms and power transformations, is the family of Box-Cox transformations with a parameter \(\lambda \ge 0\) : \[ w_t = \begin{cases} \begin{split} & \log(x_t) \quad \text{ if }\lambda=0 \\ & (sign(x_t)|x_t|^\lambda-1)/\lambda \quad \text{ otherwise } \end{split} \end{cases} \]
  • A value of \(\lambda=0\) corresponds to the multiplicative decomposition
  • A value of \(\lambda=1\) corresponds to the additive decomposition

stl() function

• For our beer production example, the following command works:

### library
library(dplyr)
library(astsa)
library(forecast)

### data
beer <- scan("data//beerprod.dat")
beer <- ts(beer, freq = 4)

### stl 1
stl.beer1 <- stl(beer, "periodic")

# plot
plot(stl.beer1)

# trend & seasonal & random
stl.beer1$time.series

       seasonal    trend   remainder
 1 Q1   8.06289 267.3569   8.9802489
 1 Q2 -41.58529 261.8710  -7.4856917
 1 Q3 -24.68456 257.1444  -5.5598227
 1 Q4  58.20698 253.8595  -3.6664533
 2 Q1   8.06289 257.4133  -3.4761521
 2 Q2 -41.58529 260.6083   8.8769848
 2 Q3 -24.68456 262.9967  -2.2121517
 2 Q4  58.20698 264.2030  -2.0099463
 3 Q1   8.06289 265.5992  -1.7621096
 3 Q2 -41.58529 265.2844   9.1008543
 3 Q3 -24.68456 262.6567  -0.9721191
 3 Q4  58.20698 259.6218  -4.4287360
 4 Q1   8.06289 260.6322  -7.2950710
 4 Q2 -41.58529 263.5667   4.8185798
 4 Q3 -24.68456 266.4737   8.1108205
 4 Q4  58.20698 268.8899 -12.7968992
 5 Q1   8.06289 270.2192   7.8179495
 5 Q2 -41.58529 272.3052  -4.2199519
 5 Q3 -24.68456 271.9999  13.0846868
 5 Q4  58.20698 273.5345 -20.3415229
 6 Q1   8.06289 274.2392  12.3978832
 6 Q2 -41.58529 277.9325  -3.7472193
 6 Q3 -24.68456 279.0722   2.8123800
 6 Q4  58.20698 279.0277   1.9653099
 7 Q1   8.06289 282.2380 -11.2008871
 7 Q2 -41.58529 285.5337   5.8516025
 7 Q3 -24.68456 288.6762   5.8083664
 7 Q4  58.20698 290.1045  -2.6114918
 8 Q1   8.06289 291.4126  -5.6754543
 8 Q2 -41.58529 295.1049   1.1803633
 8 Q3 -24.68456 300.1132   2.0713696
 8 Q4  58.20698 304.4871   0.7059037
 9 Q1   8.06289 309.2504  -3.9132616
 9 Q2 -41.58529 313.5026   0.8827153
 9 Q3 -24.68456 316.5087   8.2759085
 9 Q4  58.20698 319.8058  -8.5127401
10 Q1   8.06289 322.4975   0.2395676
10 Q2 -41.58529 326.0288   3.3565273
10 Q3 -24.68456 328.4569   2.1276459
10 Q4  58.20698 328.7079  -0.8148974
11 Q1   8.06289 328.8886  -1.7514832
11 Q2 -41.58529 331.1114  -1.5260929
11 Q3 -24.68456 335.5519  -2.5673557
11 Q4  58.20698 341.3185   2.7745080
12 Q1   8.06289 347.2008  -2.4637161
12 Q2 -41.58529 349.8189   7.8664098
12 Q3 -24.68456 353.9466  -4.3620815
12 Q4  58.20698 359.6907 -13.0976594
13 Q1   8.06289 361.9729  22.9641718
13 Q2 -41.58529 365.4701  -4.9848011
13 Q3 -24.68456 367.9136 -16.2290100
13 Q4  58.20698 370.0087  14.0842803
14 Q1   8.06289 375.6530  -0.6158554
14 Q2 -41.58529 379.5102  -6.3249053
14 Q3 -24.68456 381.3293   4.7552906
14 Q4  58.20698 384.5914   3.1015678
15 Q1   8.06289 388.8069 -10.2697510
15 Q2 -41.58529 393.2770   5.5082520
15 Q3 -24.68456 399.2517  -0.9671132
15 Q4  58.20698 403.2668   4.7261825
16 Q1   8.06289 405.5667  -4.0295442
16 Q2 -41.58529 408.1858   3.1995104
16 Q3 -24.68456 412.5123  -9.2277559
16 Q4  58.20698 414.9332  13.8597888
17 Q1   8.06289 417.3459  -6.2088201
17 Q2 -41.58529 420.2610  -1.9757578
17 Q3 -24.68456 428.1381 -10.6535322
17 Q4  58.20698 435.8692  12.0238431
18 Q1   8.06289 441.0740   9.2631448
18 Q2 -41.58529 447.1144 -18.1291496
18 Q3 -24.68456 453.0243  -1.4397012
18 Q4  58.20698 459.3832   7.4098529

# seasonally adjusted values
seasadj(stl.beer1)

       Qtr1     Qtr2     Qtr3     Qtr4
1  276.3371 254.3853 251.5846 250.1930
2  253.9371 269.4853 260.7846 262.1930
3  263.8371 274.3853 261.6846 255.1930
4  253.3371 268.3853 274.5846 256.0930
5  278.0371 268.0853 285.0846 253.1930
6  286.6371 274.1853 281.8846 280.9930
7  271.0371 291.3853 294.4846 287.4930
8  285.7371 296.2853 302.1846 305.1930
9  305.3371 314.3853 324.7846 311.2930
10 322.7371 329.3853 330.5846 327.8930
11 327.1371 329.5853 332.9846 344.0930
12 344.7371 357.6853 349.5846 346.5930
13 384.9371 360.4853 351.6846 384.0930
14 375.0371 373.1853 386.0846 387.6930
15 378.5371 398.7853 398.2846 407.9930
16 401.5371 411.3853 403.2846 428.7930
17 411.1371 418.2853 417.4846 447.8930
18 450.3371 428.9853 451.5846 466.7930

beer - stl.beer1$time.series[,1]

       Qtr1     Qtr2     Qtr3     Qtr4
1  276.3371 254.3853 251.5846 250.1930
2  253.9371 269.4853 260.7846 262.1930
3  263.8371 274.3853 261.6846 255.1930
4  253.3371 268.3853 274.5846 256.0930
5  278.0371 268.0853 285.0846 253.1930
6  286.6371 274.1853 281.8846 280.9930
7  271.0371 291.3853 294.4846 287.4930
8  285.7371 296.2853 302.1846 305.1930
9  305.3371 314.3853 324.7846 311.2930
10 322.7371 329.3853 330.5846 327.8930
11 327.1371 329.5853 332.9846 344.0930
12 344.7371 357.6853 349.5846 346.5930
13 384.9371 360.4853 351.6846 384.0930
14 375.0371 373.1853 386.0846 387.6930
15 378.5371 398.7853 398.2846 407.9930
16 401.5371 411.3853 403.2846 428.7930
17 411.1371 418.2853 417.4846 447.8930
18 450.3371 428.9853 451.5846 466.7930

stl.beer1$time.series[,-1] %>% apply(1,sum)

 [1] 276.3371 254.3853 251.5846 250.1930 253.9371 269.4853 260.7846
 [8] 262.1930 263.8371 274.3853 261.6846 255.1930 253.3371 268.3853
[15] 274.5846 256.0930 278.0371 268.0853 285.0846 253.1930 286.6371
[22] 274.1853 281.8846 280.9930 271.0371 291.3853 294.4846 287.4930
[29] 285.7371 296.2853 302.1846 305.1930 305.3371 314.3853 324.7846
[36] 311.2930 322.7371 329.3853 330.5846 327.8930 327.1371 329.5853
[43] 332.9846 344.0930 344.7371 357.6853 349.5846 346.5930 384.9371
[50] 360.4853 351.6846 384.0930 375.0371 373.1853 386.0846 387.6930
[57] 378.5371 398.7853 398.2846 407.9930 401.5371 411.3853 403.2846
[64] 428.7930 411.1371 418.2853 417.4846 447.8930 450.3371 428.9853
[71] 451.5846 466.7930

• We may conisder three arguments :
  1. s.window
     • The user must specify s.window as there is no default
     • s.window = "periodic means that the seasonal component is periodic (i.e., identical across years)
     • we may set s.window as an odd number (it should be odd) : it will be number of consecutive years to be used in estimating each value in the seasonal component. That is, the seasonal component is allowed to change over time
  2. t.window
     • Specifying t.window is optional, and a default value will be used if it is omitted.
     • we may set t.window as an odd number (it should be odd) : it will be the number of consecutive observations to be used when estimating the trend component
  3. robust
     • robust = F is the default setting.
     • robust = T will give a robust estimation for the trend and seasonal component.

### stl 2
stl.beer2 <- stl(beer, t.window=5, s.window=5, robust=T)

# plot
plot(stl.beer2)

# trend & seasonal & random
stl.beer2$time.series[,1]

         Qtr1       Qtr2       Qtr3       Qtr4
1    6.129476 -31.655616 -26.448979  51.924341
2    6.187253 -31.604769 -26.265588  51.755440
3    5.987362 -31.754047 -26.418262  50.046702
4    9.169582 -33.078276 -25.998917  48.877802
5   11.927482 -35.964547 -24.620839  51.764272
6   12.491247 -43.276144 -22.848418  55.827656
7    8.497663 -41.664406 -21.612635  57.289364
8    3.932941 -40.787327 -21.476997  58.226819
9    4.492452 -40.840407 -22.334687  58.542510
10   5.545564 -41.441322 -23.567205  59.417830
11   6.432712 -41.594603 -25.366731  60.145894
12   6.718950 -40.735446 -26.497008  61.846406
13   3.848321 -39.039060 -25.786495  62.013785
14   1.107160 -37.698541 -25.698301  61.737225
15   1.936468 -37.001483 -26.854931  61.614043
16   2.768080 -37.188031 -28.285833  63.276082
17   2.606293 -37.780802 -28.735462  64.260088
18   2.493020 -37.972076 -28.908200  64.462103

---

Section 6.2 Smoothing

• Smoothing is usually done to help us better see patterns (for example, trends) in time series.
• Generally smooth out the irregular roughness to see a clearer signal.
• For seasonal data, we might smooth out the seasonality so that we can identify the trend.
• Smoothing doesn’t provide us with a model, but it can be a good first step in describing various components of the series.
• The term filter is sometimes used to describe a smoothing procedure.
  • For instance, if the smoothed value for a particular time is calculated as a linear combination of observations for surrounding times, it might be said that we’ve applied a linear filter to the data.
• Filtering, smoothing, and forecasting can be strictly defined as follows: We want to produce estimators for \(x_t\) given the data \(\{ x_1, \dots , x_s \}\)
  • Forecasting (or prediction) : when \(t > s\)
  • Filtering : when \(t = s\)
  • Smoothing : when \(t < s\)

---

6.2.1 Moving averages

• A moving average method is one of the simplest smoothing methods
• The classical method of time series decomposition uses a moving average method to estimate the trend component.
• To take away seasonality from a series so we can better see trend, we would use a moving average with a length = seasonal span.
  • Thus in the smoothed series, each smoothed value has been averaged across all seasons.
  • For this, there is 2 approaches :
    1. “one-sided” moving average : you average all values before the current time
    2. “two-sided” (or centered) moving average : you average values both before and after the current time.
• (Example) “one-sided” moving average
  • For a quarterly time-series data (Frequency = 4), we could define a smoothed value for time \(t\) as \((x_t + x_{t-1} + x_{t-2} + x_{t-3})/4\), the average of this time and the previous 3 quarters.
  • In R code, this will be a one-sided filter.
• “two-sided” (or centered) moving average :
  • Moving average of order \(m=2k+1\) : \[ \hat{T}_t = {1 \over m} \sum_{=-k}^{k}{x_{t+j}} \]
    • (Example) For a daily time-series data with a weekly seasonality (Frequency = 7), moving average of order \(m=7\) can be used : \[ \hat{T}_t = {1 \over 7}(x_{t-3} + x_{t-2} + \dots + x_{t+3}) \]
  • Moving average of order \(m=2k\) : We apply a moving average to a moving average \[ \hat{T}_t = {1 \over 2m} x_{t-k} + {1 \over m} \sum_{-k+1}^{k-1}{x_{t+j}} + {1 \over 2m} x_{t+k} \]
    • (Example) For a quarterly time-series data (Frequency = 4), we have the following :
      • Here’s how the centered moving average for \(t = 3\) would be calculated.
      • Average the observed data values at times 1 to 4: \[ {1\over 4}(x_1 + x_2 + x_3 + x_4) \]
      • Then average those two averages: \[\begin{equation} \begin{split} & {1\over 2} \left( {1\over 4}(x_1 + x_2 + x_3 + x_4) + {1\over 4}(x_2 + x_3 + x_4 + x_5) \right) \\ & = {1 \over 8}x_1 + {1 \over 4}x_2 + {1 \over 4}x_3 + {1 \over 4}x_4 + {1 \over 8}x_5 \end{split} \end{equation}\]
    • (Exmple; Cont.) To smooth away seasonality in quarterly data, in order to identify trend, the usual convention is to use the moving average smoothed at time \(t\) is \[ {1 \over 8}x_{t-2} + {1 \over 4}x_{t-1} + {1 \over 4}x_t + {1 \over 4}x_{t+1} + {1 \over 8}x_{t+2} \]
    • (Example) To smooth away seasonality in monthly data, in order to identify trend, the usual convention is to use the moving average smoothed at time \(t\) is \[ {1 \over 24}x_{t-6} + {1 \over 12}\sum_{j=t-5}^{t+5}{x_{j}} + {1 \over 24}x_{t+6} \]
• For the filter command in R, we’ll specify a two-sided filter when we want to use values that come both before and after the time for which we’re smoothing.

---

Example 3

• We look at a series of quarterly beer production in Australia.
• The following R code creates a smoothed series that lets us see the trend pattern, and plots this trend pattern on the same graph as the time series.
• The second command creates and stores the smoothed series in the object called trendpattern.
• Note that Within the filter command, the parameter named filter gives the coefficients for our smoothing and sides = 2 causes a centered smooth to be calculated.

beer <- scan("data//beerprod.dat")
trendpattern <- stats::filter(beer, filter = c(1/8, 1/4, 1/4, 1/4, 1/8), sides=2)
plot(beer, type="b", main = "moving average annual trend")
lines(trendpattern, col=4, lwd=2)

• We might subtract the trend pattern from the data values to get a better look at seasonality. Here’s how that would be done:

seasonals <- beer - trendpattern
plot(seasonals, type = "b", main = "Seasonal pattern for beer production")

• Another possibility for smoothing series to see trend is the one-sided filter
  • With this, the smoothed value is the average of the past year.

trendpattern2 <- stats::filter(beer, filter = c(1/4, 1/4, 1/4, 1/4), sides=1)

# plot 1
plot(beer, type="b", main = "moving average annual trend2")
lines(trendpattern2, col=4, lwd=2)

# seasonality
seasonals2 <- beer - trendpattern2
plot(seasonals2, type = "b", main = "Seasonal pattern2 for beer production")

---

Example 4

• We look at a monthly series of U.S. Unemployment for 1948-1978. Here’s a smoothing done to look at the trend.

trendunemploy <- stats::filter(unemp, filter = c(1/24,1/12,1/12,1/12,1/12,1/12,1/12,1/12,1/12,1/12,1/12,1/12,1/24), sides = 2)
trendunemploy <- ts(trendunemploy, start = c(1948,1), freq = 12)
plot(unemp, main="Trend in U.S. Unemployment, 1948-1978", xlab = "Year")
lines(trendunemploy, col=4, lwd=2)

---

Non-Seasonal Series

• For non-seasonal series, you aren’t bound to smooth over any particular span.
• For smoothing, you should experiment with moving averages of different spans.
  • Those spans of time could be relatively short.
  • The objective is to knock off the rough edges to see what trend or pattern might be there.

---

Other Smoothing Methods

• Section 2.3 of the optional textbook describes several sophisticated and useful alternatives to moving average smoothing.
• Here, we conisder the method of lowess (locally weighted regression), which may be the most widely used one.
• The following plot is the smoothed trend line for the U.S. Unemployment series.

unemp <- ts(unemp, start = c(1948,1), freq=12)

# lowess 1
unemp_lowess <- lowess(unemp, f = 2/3)
plot(unemp, main ="Lowess smoothing 1 of U.S. Unemployment Trend")
lines(unemp_lowess, lwd=2, col=4)

# lowess 2
unemp_lowess <- lowess(unemp, f = 0.1)
plot(unemp, main ="Lowess smoothing 2 of U.S. Unemployment Trend")
lines(unemp_lowess, lwd=2, col=4)

# lowess 2
unemp_lowess <- lowess(unemp, f = 1)
plot(unemp, main ="Lowess smoothing 3 of U.S. Unemployment Trend")
lines(unemp_lowess, lwd=2, col=4)

---

6.2.2 Single Exponential Smoothing

• The simplest one of the exponentially smoothing methods is the single (or simple) exponential smoothing (SES).

• This method is suitable for forecasting data with no clear trend or seasonal pattern.

  • Although the method is called a smoothing method, it’s principally used for short run forecasting.

• For forecasting, we may consider the following naive method, where all forecasts for the future are equal to the last observed value of the series : \[ \hat{x}_{t+h|t} = x_t \] for \(h=1,2,\dots\)

  • Here, the naive method assumes that the most recent observation is the only important one, and all previous observations provide no information for the future.
  • This can be thought of as a weighted average where all of the weight is given to the last observation.

• We may consider another method as follows, where all future forecasts are equal to a simple average of the observed data : \[ \hat{x}_{t+h|T} = {1\over t} \sum_{s=1}^{t}{x_s} \] for \(h=1,2,\dots\)

  • Hence, the average method assumes that all observations are of equal importance, and gives them equal weights when generating forecasts.

• The single exponential smoothing method is something between these two extremes above : \[ \hat{x}_{t+1} = \alpha x_t + (1-\alpha) \hat{x}_t \]

  • We forecast the value of \(x\) at time \(t+1\) to be a weighted combination of the observed value at time \(t\) and the forecasted value at time \(t\).
  • Actually, this is a weighted average method, where the weights decrease exponentially as observations come from further in the past (Why?) : \[ \hat{x}_{t+1} = \alpha x_t + \alpha(1-\alpha) x_{t-1} + \alpha(1-\alpha)^2 x_{t-2} + \dots \]

• The value of \(0 \le \alpha \le 1\) is called the smoothing constant.

  • \(\alpha= 0.2\) is a popular default choice of programs.
    • This puts a weight of \(0.2\) on the most recent observation and a weight of \(1 − 0.2 = 0.8\) on the most recent forecast.
  • With a relatively small value of \(\alpha\), the smoothing will be relatively more extensive.
  • With a relatively large value of \(\alpha\), the smoothing is relatively less extensive as more weight will be put on the observed value.

• This is simple one-step ahead forecasting method that at first glance seems not to require a model for the data.

• In fact, this method is equivalent to the use of an ARIMA\((0,1,1)\) model with no constant. (Why?)

  • The optimal procedure is to fit an ARIMA\((0,1,1)\) model to the observed dataset and use the results to determine the value of \(\alpha\).
  • This is “optimal” in the sense of creating the best \(\alpha\) for the data already observed.

• Although the goal is smoothing and one step ahead forecasting, the equivalence to the ARIMA\((0,1,1)\) model does bring up a good point.

  • We shouldn’t blindly apply exponential smoothing because the underlying process might not be well modeled by an ARIMA\((0,1,1)\).

---

ARIMA\((0,1,1)\) and Exponential Smoothing Equivalence

• Consider an ARIMA\((0,1,1)\) with mean \(\mu=0\) for the first differences, \(x_t - x_{t-1}\):
• The model is \(x_t-x_{t-1}=w_t + \theta_1w_{t-1}\).
  • Equivalently, \(x_t = x_{t-1} + w_t + \theta_1w_{t-1}\).
• To forecast at time \(t+1\), we consider \(x_{t+1} = x_t + w_{t+1} + \theta_1 w_t\).
  • Because \(w_{t+1} = x_{t+1} - \hat{x}_{t+1}\), we have \[\begin{equation} \begin{split} \hat{x}_{t+1} & = x_t + \theta_1 w_t \\ & = x_t + \theta_1 (x_t - \hat{x}_t) \\ & = (1 + \theta_1)x_t - \theta_1 \hat{x}_t \\ \end{split} \end{equation}\]
• If we let \(\alpha = (1+ \theta_1)\) and thus we have \(-(\theta_1) = 1−\alpha\)

---

Why the Method is Called Exponential Smoothing

• Starting with \(\hat{x}_{t+1} = \alpha x_t + (1-\alpha) \hat{x}_{t}\), we have the following : \[\begin{equation} \begin{split} \hat{x}_{t+1} & = \alpha x_t + (1-\alpha) \hat{x}_{t} \\ \hat{x}_{t} & = \alpha x_{t-1} + (1-\alpha) \hat{x}_{t-1} \\ \vdots & \\ \hat{x}_{4} & = \alpha x_{3} + (1-\alpha) \hat{x}_{3} \\ \hat{x}_{3} & = \alpha x_{2} + (1-\alpha) \hat{x}_{2} \\ \hat{x}_{2} & = \alpha x_{1} + (1-\alpha) \hat{x}_{1} \\ \end{split} \end{equation}\]

• For \(\hat{x}_{t+1} = \alpha x_t + (1-\alpha) \hat{x}_{t}\), we can substitute for \(\hat{x}_t\) : \[\begin{equation} \begin{split} \hat{x}_{t+1} & = \alpha x_t + (1-\alpha)[\alpha x_{t-1} + (1-\alpha)\hat{x}_{t-1}] \\ & = \alpha x_t + \alpha(1-\alpha)x_{t-1} + (1-\alpha)^2 \hat{x}_{t-1} \\ \end{split} \end{equation}\]

• Continue in this fashion by successively substituting for the forecasted value on the right side of the equation. This leads to: \[\begin{equation} \begin{split} \hat{x}_{t+1} = & \alpha x_t + \alpha(1-\alpha)x_{t-1} + \alpha(1-\alpha)^2 x_{t-2} + \dots + \\ & \alpha(1-\alpha)^j x_{t-j} + \dots + \alpha(1-\alpha)^{t-1} x_{1} + (1-\alpha)^{t} \hat{x}_{1} \\ & \sum_{j=0}^{t-1}{ \alpha (1-\alpha)^j x_{t-j} }+ (1-\alpha)^{t} \hat{x}_{1} \\ \end{split} \end{equation}\]

• The above equation shows that the forecasted value is a weighted average of all past values of the series, with exponentially changing weights as we move back in the series.

---

Optimal Exponential Smoothing

• Note that a single xxponential smoothing require the smoothing parameter \(\alpha\) and the initial value \(x_0 = \hat{x}_{1}\) to be chosen.

  • Note that the selection of \(x_0\) has less effect on the single exponential smoothing as \(t\) gets bigger.

• There is no formally correct procedure for choosing the smoothing parameter \(\alpha\) and the initial value \(x_0\).

• One approach for choosing appropriate smoothing parameter \(\alpha\) and the initial value \(x_0\) : we find we find the values of the unknown parameter and the initial value that minimize \[ \text{SSE} = \sum_{s=1}^{t}{(x_t - \hat{x}_t)^2} = \sum_{s=1}^{t}{e_t^2} \]

• Or we may consider this approach based the equivalence with an ARIMA\((0,1,1)\) :

  • Basically, we just fit an ARIMA\((0,1,1)\) to the data and determine the \(\alpha\) coefficient.
  • We can examine the fit of the smooth by comparing the predicted values to the actual series.

---

Components form

• An alternative representation is the component form: The component form of simple exponential smoothing is given by: \[\begin{equation} \begin{split} \text{Forecasting equation : } & \hat{x}_{t+h|t} & = l_t \\ \text{Smoothing equation : } & l_t & = \alpha x_t + (1-\alpha)l_{t-1} \\ \end{split} \end{equation}\]
  • where \(l_t\) is an estimate of the level (or the smoothed value) of the series at time \(t\).
• Note that a simple exponential smoothing includes only the level component \(l_t\).
  • Other advanced methods may also include a trend \(b_t\) and a seasonal component \(s_t\).

---

Flat forecasts

• Simple exponential smoothing has a “flat” forecast function : \[ \hat{x}_{t+h|t} = \hat{x}_{t+1|t} = \hat{x}_{t+1} = l_t \] for \(h=2,3\dots\)
  • That is, all forecasts take the same value, equal to the last level component
  • Therefore, these forecasts will only be suitable if the time-series has no trend or seasonal component.

---

Example 5

• We consider oil price index data:

### library
library(dplyr)
library(astsa)
library(forecast)
library(fpp2)

### data
plot(oil, main="oil production")

### ses
result <- ses(oil, h=5)
result

     Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
2014       542.3412 479.4802 605.2022 446.2037 638.4788
2015       542.3412 453.4468 631.2356 406.3890 678.2935
2016       542.3412 433.4701 651.2124 375.8372 708.8453
2017       542.3412 416.6287 668.0537 350.0805 734.6019
2018       542.3412 401.7911 682.8914 327.3883 757.2941

plot(result)

### ARIMA(0,1,1)
fit0 <- sarima(oil, 0,1,1, no.constant=T)

initial  value 3.882308 
iter   2 value 3.866729
iter   3 value 3.866238
iter   4 value 3.866237
iter   4 value 3.866237
final  value 3.866237 
converged
initial  value 3.866482 
iter   2 value 3.866478
iter   3 value 3.866478
iter   3 value 3.866478
iter   3 value 3.866478
final  value 3.866478 
converged

fit0 # 0.1627

$fit

Call:
arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = S), 
    include.mean = !no.constant, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc, 
        REPORT = 1, reltol = tol))

Coefficients:
         ma1
      0.1627
s.e.  0.1310

sigma^2 estimated as 2281:  log likelihood = -253.7,  aic = 511.4

$degrees_of_freedom
[1] 47

$ttable
    Estimate    SE t.value p.value
ma1   0.1627 0.131  1.2414  0.2206

$AIC
[1] 10.65417

$AICc
[1] 10.65598

$BIC
[1] 10.73213

---

6.2.3 Double Exponential Smoothing (Holt’s linear trend method)

• Double exponential smoothing (Holt’s linear trend method) might be used when there’s trend (either long run or short run), but no seasonality.
• The Holt’s linear trend method can be expressed as follows : \[\begin{equation} \begin{split} \text{Forecasting equation : } & \hat{x}_{t+h|t} & = l_t \\ \text{Level equation : } & l_t & = \alpha x_t + (1-\alpha)(l_{t-1}+b_{t-1}) \\ \text{Trend equation : } & b_t & = \beta^* (l_t-l_{t-1}) + (1-\beta^*)b_{t-1} \\ \end{split} \end{equation}\]
  • where \(l_t\) is an estimate of the level (or the smoothed value) of the series at time \(t\),
  • \(b_t\) denotes an estimate of the trend (slope) of the series at time \(t\)
  • \(0 \le \alpha \le 1\)is the smoothing parameter for the level
  • \(0 \le \beta^* \le 1\)is the smoothing parameter for the trend

### library
library(dplyr)
library(astsa)
library(forecast)
library(fpp2)

### data
plot(oil, main="oil production")

### holt
result <- holt(oil, h=5)
result

     Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
2014       554.0762 489.1032 619.0491 454.7086 653.4438
2015       565.8091 473.9255 657.6927 425.2852 706.3330
2016       577.5420 465.0051 690.0789 405.4317 749.6524
2017       589.2750 459.3236 719.2263 390.5315 788.0184
2018       601.0079 455.7116 746.3041 378.7964 823.2193

plot(result)

---

6.2.4 ets() function in R

• There are also further extensions for exponential smoothing methods, including the seasonality component and nonlinear trend and so on.
• In R, we can perform the various exponential smoothing methods and model selection by using the est() function in the forecast package.
  • By default it uses the AICc to select an appropriate model, although other information criteria can be selected.
• (Example 1)

### library
library(dplyr)
library(astsa)
library(forecast)
library(fpp2)

### data
plot(oil, main="oil production")

### ets
fit <- ets(oil)

### forecasting
fit %>% forecast(h=10) %>% plot()

• (Example 2)

### library
library(dplyr)
library(astsa)
library(forecast)
library(fpp2)

### data
plot(austourists, main="International tourist visitor nights in Australia")

### ets
fit <- ets(austourists)

### forecasting
fit %>% forecast(h=10) %>% plot()

• ETS(\(\cdot,\cdot,\cdot\)) indicates (Error, Trend, Seasonal) where
  • Error = \(\{ A, M \}\)
    • Additive error / Multiplicative error
  • Trend = \(\{ N, A, A_d \}\)
    • No trend / Linear trend / Damped trend
  • Seasonal = \(\{ N, A, M \}\)
    • No Seasonality / Additive Seasonality / Multiplicative Seasonality