Base R

Alexander Häußer

Source: vignettes/vignette_01_baseR.Rmd

Load package 

library(echos)
#> Registered S3 method overwritten by 'quantmod':
#>   method            from
#>   as.zoo.data.frame zoo

Prepare dataset

In a first example, we want to model the well-known AirPassenger time series (ts object). The dataset contains monthly totals of international airline passengers (in thousands) from January 1949 to December 1960 with 144 observations in total. The first 132 observations are used for model training (n_train) and last 12 observations are used for testing (n_ahead). xtrain and xtest are numeric vectors containing the training and testing data.

# Forecast horizon
n_ahead <- 12
# Number of observations (total)
n_obs <- length(AirPassengers)
# Number of observations (training data)
n_train <- n_obs - n_ahead

# Prepare train and test data as numeric vectors
xtrain <- AirPassengers[(1:n_train)]
xtest <- AirPassengers[((n_train+1):n_obs)]

xtrain
#>   [1] 112 118 132 129 121 135 148 148 136 119 104 118 115 126 141 135 125 149
#>  [19] 170 170 158 133 114 140 145 150 178 163 172 178 199 199 184 162 146 166
#>  [37] 171 180 193 181 183 218 230 242 209 191 172 194 196 196 236 235 229 243
#>  [55] 264 272 237 211 180 201 204 188 235 227 234 264 302 293 259 229 203 229
#>  [73] 242 233 267 269 270 315 364 347 312 274 237 278 284 277 317 313 318 374
#>  [91] 413 405 355 306 271 306 315 301 356 348 355 422 465 467 404 347 305 336
#> [109] 340 318 362 348 363 435 491 505 404 359 310 337 360 342 406 396 420 472
#> [127] 548 559 463 407 362 405
xtest
#>  [1] 417 391 419 461 472 535 622 606 508 461 390 432

Train ESN model

The function train_esn() is used to automatically train an ESN to the input data xtrain, where the output xmodel is an object of class esn. The object xmodel is a list containing the actual and fitted values, residuals, the internal states states_train, estimated coefficients from the ridge regression estimation, hyperparameters, etc. We can summarize the model by using the generic S3 method summary() to get detailed information on the trained model.

# Train ESN model
xmodel <- train_esn(y = xtrain)

# Plot actual and fitted values
plot(xmodel$actual, type = "l")
lines(xmodel$fitted, col = "steelblue", lwd = 2)


# Summarize model
summary(xmodel)
#> 
#> --- Layers ----------------------------------------------------- 
#>  n_inputs  =  1 
#>  n_states  =  52 
#>  n_outputs =  1 
#> 
#> --- Meta --------------------------------------------------- 
#>  lags     =  1 
#>  n_diff   =  1 
#>  n_models =  104 
#> 
#> --- Scaling ----------------------------------------------------
#> scale_inputs = [-0.5, 0.5]
#> scale_win    = [-0.5, 0.5]
#> scale_wres   = [-0.5, 0.5]
#> 
#> --- Hyperparameters -------------------------------------------- 
#>  alpha   =  1 
#>  rho     =  1 
#>  density =  0.5

From the output above, we get the following information about the trained model:

  • Layers
    • n_inputs: Number of input variables
    • n_states: Number of internal states (i.e., the reservoir size). The internal states are the actual predictor variables in the model.
    • n_outputs: Number of output variables (i.e., response variables)
  • Meta
    • lags: Lags of the output variable, which are used as model input
    • n_diff: Number of differences required to achieve (weak-) stationarity of the input training data
    • n_models: Number of models evaluated during the random search optimization to find the regularization parameter lambda
  • Scaling
    • scale_inputs: Input training data are scaled to the interval (-0.5, 0.5)
    • scale_win: Input weights matrix is drawn from a random uniform distribution with interval (-0.5, 0.5)
    • scale_wres: Reservoir weight matrix is drawn from a random uniform distribution with interval (-0.5, 0.5)
  • Hyperparameters:
    • alpha: Leakage rate (smoothing parameter)
    • rho: Spectral radius for scaling the reservoir weight matrix
    • lambda: Regularization parameter for the ridge regression estimation
    • density: The density of the reservoir weight matrix

Forecast ESN model

The function forecast_esn() is used to forecast the trained model xmodel for n_ahead steps into the future. The output xfcst is a list of class forecast_esn containing the point forecasts, actual and fitted values, the forecast horizon n_ahead and the model specification model_spec. We can use the generic S3 method plot() to visualize the point forecast within xfcst and add the holdout test data xtest.

# Forecast ESN model
xfcst <- forecast_esn(xmodel, n_ahead = n_ahead)
xfcst
#> $point
#>  [1] 417.1882 394.6582 442.8875 432.9793 466.5263 522.0339 594.8509 601.8089
#>  [9] 501.4787 445.6808 400.5648 438.0486
#> 
#> $actual
#>   [1]  NA  NA  NA  NA  NA  NA  NA  NA 136 119 104 118 115 126 141 135 125 149
#>  [19] 170 170 158 133 114 140 145 150 178 163 172 178 199 199 184 162 146 166
#>  [37] 171 180 193 181 183 218 230 242 209 191 172 194 196 196 236 235 229 243
#>  [55] 264 272 237 211 180 201 204 188 235 227 234 264 302 293 259 229 203 229
#>  [73] 242 233 267 269 270 315 364 347 312 274 237 278 284 277 317 313 318 374
#>  [91] 413 405 355 306 271 306 315 301 356 348 355 422 465 467 404 347 305 336
#> [109] 340 318 362 348 363 435 491 505 404 359 310 337 360 342 406 396 420 472
#> [127] 548 559 463 407 362 405
#> 
#> $fitted
#>   [1]        NA        NA        NA        NA        NA        NA        NA
#>   [8]        NA 133.71643 110.30808  97.99022 121.00545 115.14693 126.53016
#>  [15] 157.52189 132.94333 119.78278 161.34975 193.59773 173.76932 143.85410
#>  [22] 107.56646 105.13810 155.56655 149.80988 160.86024 198.54238 147.81409
#>  [29] 164.43319 197.62710 225.94828 197.05461 167.66886 141.42018 124.95530
#>  [36] 195.90012 173.94581 188.78576 221.17353 172.36502 185.81144 233.87344
#>  [43] 253.80249 242.09268 185.71449 162.28386 157.26375 217.82387 205.98666
#>  [50] 197.74646 261.03987 215.40821 240.28469 268.54261 279.86910 276.93788
#>  [57] 212.74332 184.36646 163.91128 226.25235 215.92186 191.58228 289.17804
#>  [64] 226.18730 244.44391 283.04411 328.54638 306.08284 220.99874 194.32638
#>  [71] 174.07617 260.50090 239.72264 219.51924 309.30297 264.75397 270.50592
#>  [78] 350.00199 399.61293 335.73425 262.16156 235.00641 207.03053 310.60081
#>  [85] 287.16356 272.65522 350.18168 309.04424 318.45190 415.49340 464.43369
#>  [92] 392.70558 304.88560 261.77706 236.78694 342.73360 327.14560 292.43738
#>  [99] 397.08165 344.62531 363.33019 473.10273 505.71933 457.04417 345.84712
#> [106] 291.90008 266.42954 369.29088 353.07476 303.67445 418.08114 343.69856
#> [113] 382.50629 507.63837 550.88035 507.75092 337.10651 301.79449 268.94196
#> [120] 370.79001 374.94003 310.78017 458.56659 377.22597 441.29224 532.46027
#> [127] 607.88392 565.09347 369.72182 353.55294 310.97893 429.61063
#> 
#> $n_ahead
#> [1] 12
#> 
#> $model_spec
#> [1] "ESN(1,52,1)"
#> 
#> attr(,"class")
#> [1] "forecast_esn"

# Plot forecast and test data
plot(xfcst, test = xtest)