These are the different algorithms that can used in modeling multiple channel EEG time series for neuroscience applications. The data set is from the description at www.neuronalarchitects.com/NeuralCipher under the EEG Data sets for the Epilepsy center in Born, Germany. It is a work in progress.
Step 1. Load libraries.
library(tsDyn)
## Loading required package: mgcv
## This is mgcv 1.7-24. For overview type 'help("mgcv-package")'.
## Loading required package: Matrix
## Loading required package: lattice
## Loading required package: snow
## Loading required package: mnormt
## Loading required package: foreach
## Loading required package: MASS
## Loading required package: nlme
library(tseriesChaos)
## Loading required package: deSolve
library(BAYSTAR)
## Error: there is no package called 'BAYSTAR'
library(coda)
library(Metrics)
## Attaching package: 'Metrics'
## The following object is masked from 'package:tsDyn':
##
## mse
Step 2. Load the Sample Data and choose a variable.
setwd("C:/CromwellWorkshp/")
## Error: cannot change working directory
EEG_Data <- read.csv("EEG_test3.csv", head = TRUE)
attach(EEG_Data)
Step 3. Data Visualization of each of the channels.
plot(EEG_Data)
plot(F1, type = "l")
plot(F2, type = "l")
plot(F3, type = "l")
plot(F4, type = "l")
plot(F5, type = "l")
plot(F6, type = "l")
plot(F7, type = "l")
plot(F8, type = "l")
plot(F9, type = "l")
plot(F10, type = "l")
Step 4. Translations and the selection of F3
# Translate to the F3 designation.
tsF <- ts(F3)
# Remove observations from the beginning
tsF <- window(tsF, start = 10, end = 4100)
## Warning: 'end' value not changed
Step 5. Specification and Estimation of the Models.
mod <- list()
# SETAR
grid <- selectSETAR(tsF, m = 3, thDelay = 1, trim = 0.15, criterion = "AIC")
## Using maximum autoregressive order for low regime: mL = 3
## Using maximum autoregressive order for high regime: mH = 3
## Searching on 154 possible threshold values within regimes with sufficient ( 15% ) number of observations
## Searching on 1386 combinations of thresholds ( 154 ), thDelay ( 1 ), mL ( 3 ) and MM ( 3 )
##
## 1 T: Trim not respected: 0.8536 0.1464 from th: 52
## 1 T: Trim not respected: 0.8536 0.1464 from th: 52
## 1 T: Trim not respected: 0.8536 0.1464 from th: 52
## 1 T: Trim not respected: 0.8536 0.1464 from th: 52
## 1 T: Trim not respected: 0.8536 0.1464 from th: 52
## 1 T: Trim not respected: 0.8536 0.1464 from th: 52
## 1 T: Trim not respected: 0.8536 0.1464 from th: 52
## 1 T: Trim not respected: 0.8536 0.1464 from th: 52
## 1 T: Trim not respected: 0.8536 0.1464 from th: 52
print(grid)
## Results of the grid search for 1 threshold
## thDelay mL mH th AIC
## 1 1 3 3 47 17207
## 2 1 3 3 51 17207
## 3 1 3 3 40 17208
## 4 1 3 3 49 17208
## 5 1 3 3 48 17208
## 6 1 3 3 52 17209
## 7 1 3 3 38 17210
## 8 1 3 3 41 17210
## 9 1 3 3 36 17211
## 10 1 3 3 50 17212
mod[["setar"]] <- setar(tsF, m = 3, mL = 1, mH = 2, thDelay = 1)
##
## 1 T: Trim not respected: 0.8536 0.1464 from th: 52
## Warning: Possible unit root in the low regime. Roots are: 0.9262
summary(mod[["setar"]])
##
## Non linear autoregressive model
##
## SETAR model ( 2 regimes)
## Coefficients:
## Low regime:
## phiL.1 const L
## 1.08 16.99
##
## High regime:
## phiH.1 phiH.2 const H
## 1.8229 -0.8637 -1.0538
##
## Threshold:
## -Variable: Z(t) = + (0) X(t)+ (1)X(t-1)+ (0)X(t-2)
## -Value: -101
## Proportion of points in low regime: 15.18% High regime: 84.82%
##
## Residuals:
## Min 1Q Median 3Q Max
## -210.96 -4.53 0.17 5.13 43.41
##
## Fit:
## residuals variance = 91.42, AIC = 18471, MAPE = 31.04%
##
## Coefficient(s):
##
## Estimate Std. Error t value Pr(>|t|)
## const L 16.9858 1.7074 9.95 < 2e-16 ***
## phiL.1 1.0797 0.0125 86.27 < 2e-16 ***
## const H -1.0538 0.1633 -6.45 1.2e-10 ***
## phiH.1 1.8229 0.0109 167.88 < 2e-16 ***
## phiH.2 -0.8637 0.0110 -78.51 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold
## Variable: Z(t) = + (0) X(t) + (1) X(t-1)+ (0) X(t-2)
##
## Value: -101
Step 6. Apply of the model selection criterion for use.
sapply(mod, AIC)
## setar
## 18471
sapply(mod, MAPE)
## setar
## 0.3104
summary(mod[["setar"]])
##
## Non linear autoregressive model
##
## SETAR model ( 2 regimes)
## Coefficients:
## Low regime:
## phiL.1 const L
## 1.08 16.99
##
## High regime:
## phiH.1 phiH.2 const H
## 1.8229 -0.8637 -1.0538
##
## Threshold:
## -Variable: Z(t) = + (0) X(t)+ (1)X(t-1)+ (0)X(t-2)
## -Value: -101
## Proportion of points in low regime: 15.18% High regime: 84.82%
##
## Residuals:
## Min 1Q Median 3Q Max
## -210.96 -4.53 0.17 5.13 43.41
##
## Fit:
## residuals variance = 91.42, AIC = 18471, MAPE = 31.04%
##
## Coefficient(s):
##
## Estimate Std. Error t value Pr(>|t|)
## const L 16.9858 1.7074 9.95 < 2e-16 ***
## phiL.1 1.0797 0.0125 86.27 < 2e-16 ***
## const H -1.0538 0.1633 -6.45 1.2e-10 ***
## phiH.1 1.8229 0.0109 167.88 < 2e-16 ***
## phiH.2 -0.8637 0.0110 -78.51 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold
## Variable: Z(t) = + (0) X(t) + (1) X(t-1)+ (0) X(t-2)
##
## Value: -101
plot(mod[["setar"]])
Step 7. Examine the Forecasting Results and Change Parameter Specifications.
set.seed(10)
mod.test <- list()
x.train <- window(tsF, end = 4000)
x.test <- window(tsF, start = 4001)
mod.test[["linear"]] <- linear(x.train, m = 2)
mod.test[["setar"]] <- setar(x.train, m = 3, mL = 1, thDelay = 1)
## Warning: Possible unit root in the low regime. Roots are: 0.9241
# mod.test[['lstar']] <- lstar(x.train, m=2, thDelay=1, trace=FALSE,
# control=list(maxit=1e5))
mod.test[["nnet"]] <- nnetTs(x.train, m = 2, size = 3, control = list(maxit = 1e+05))
## # weights: 13
## initial value 22759704.429118
## iter 10 value 13039962.058761
## iter 20 value 10403672.629763
## iter 30 value 3864258.451380
## iter 40 value 3259847.975200
## iter 50 value 2837117.644913
## iter 60 value 2014641.110852
## iter 70 value 1327360.952112
## iter 80 value 1026222.566984
## iter 90 value 641064.227814
## iter 100 value 408534.066037
## iter 110 value 383701.468720
## iter 120 value 305518.873566
## iter 130 value 303078.432032
## iter 140 value 281112.535996
## iter 150 value 267238.874522
## iter 160 value 254772.973663
## iter 170 value 253726.265338
## iter 180 value 252841.285166
## iter 190 value 252685.040968
## iter 200 value 252665.340883
## iter 210 value 252617.376535
## final value 252576.858455
## converged
mod.test[["aar"]] <- aar(x.train, m = 2)
## Warning: the matrix is either rank-deficient or indefinite
## Warning: the matrix is either rank-deficient or indefinite
## Warning: the matrix is either rank-deficient or indefinite
## Warning: the matrix is either rank-deficient or indefinite
frc.test <- lapply(mod.test, predict, n.ahead = 100)
frc.test
## $linear
## Time Series:
## Start = 4001
## End = 4100
## Frequency = 1
## [1] -130.3774 -147.4606 -156.8637 -159.4170 -156.0982 -147.9701 -136.1233
## [8] -121.6278 -105.4922 -88.6329 -71.8508 -55.8172 -41.0670 -27.9983
## [15] -16.8777 -7.8504 -0.9536 3.8684 6.7469 7.8717 7.4738
## [22] 5.8101 3.1482 -0.2458 -4.1180 -8.2352 -12.3910 -16.4100
## [29] -20.1501 -23.5031 -26.3935 -28.7767 -30.6358 -31.9777 -32.8295
## [36] -33.2339 -33.2452 -32.9249 -32.3385 -31.5519 -30.6288 -29.6282
## [43] -28.6033 -27.5995 -26.6545 -25.7974 -25.0493 -24.4234 -23.9259
## [50] -23.5568 -23.3109 -23.1787 -23.1476 -23.2029 -23.3287 -23.5088
## [57] -23.7273 -23.9692 -24.2209 -24.4706 -24.7086 -24.9268 -25.1197
## [64] -25.2833 -25.4156 -25.5161 -25.5858 -25.6267 -25.6417 -25.6342
## [71] -25.6082 -25.5675 -25.5163 -25.4581 -25.3966 -25.3347 -25.2750
## [78] -25.2196 -25.1701 -25.1275 -25.0925 -25.0654 -25.0459 -25.0338
## [85] -25.0282 -25.0284 -25.0335 -25.0425 -25.0544 -25.0683 -25.0833
## [92] -25.0986 -25.1135 -25.1275 -25.1402 -25.1512 -25.1604 -25.1677
## [99] -25.1730 -25.1766
##
## $setar
## Time Series:
## Start = 4001
## End = 4100
## Frequency = 1
## [1] -132.26 -125.89 -119.00 -111.54 -103.47 -94.74 -85.29 -76.10
## [9] -67.70 -60.31 -53.99 -48.69 -44.29 -40.69 -37.76 -35.37
## [17] -33.44 -31.88 -30.61 -29.59 -28.76 -28.08 -27.54 -27.09
## [25] -26.73 -26.43 -26.19 -25.99 -25.83 -25.70 -25.60 -25.51
## [33] -25.44 -25.38 -25.33 -25.29 -25.26 -25.23 -25.21 -25.19
## [41] -25.18 -25.17 -25.16 -25.15 -25.15 -25.14 -25.14 -25.13
## [49] -25.13 -25.13 -25.13 -25.12 -25.12 -25.12 -25.12 -25.12
## [57] -25.12 -25.12 -25.12 -25.12 -25.12 -25.12 -25.12 -25.12
## [65] -25.12 -25.12 -25.12 -25.12 -25.12 -25.12 -25.12 -25.12
## [73] -25.12 -25.12 -25.12 -25.12 -25.12 -25.12 -25.12 -25.12
## [81] -25.12 -25.12 -25.12 -25.12 -25.12 -25.12 -25.12 -25.12
## [89] -25.12 -25.12 -25.12 -25.12 -25.12 -25.12 -25.12 -25.12
## [97] -25.12 -25.12 -25.12 -25.12
##
## $nnet
## Time Series:
## Start = 4001
## End = 4100
## Frequency = 1
## [1] -131.238 -148.848 -157.698 -158.628 -153.206 -143.249 -130.439
## [8] -116.117 -101.251 -86.492 -72.255 -58.801 -46.300 -34.875
## [15] -24.624 -15.633 -7.975 -1.699 3.176 6.672 8.853
## [22] 9.823 9.720 8.704 6.947 4.621 1.888 -1.102
## [29] -4.219 -7.354 -10.415 -13.331 -16.046 -18.520 -20.728
## [36] -22.653 -24.290 -25.639 -26.710 -27.514 -28.070 -28.397
## [43] -28.519 -28.460 -28.243 -27.894 -27.437 -26.896 -26.293
## [50] -25.649 -24.985 -24.316 -23.658 -23.025 -22.428 -21.875
## [57] -21.373 -20.928 -20.541 -20.216 -19.950 -19.742 -19.591
## [64] -19.491 -19.439 -19.429 -19.456 -19.515 -19.600 -19.706
## [71] -19.827 -19.959 -20.097 -20.237 -20.376 -20.510 -20.637
## [78] -20.755 -20.863 -20.959 -21.043 -21.113 -21.172 -21.217
## [85] -21.252 -21.275 -21.288 -21.291 -21.287 -21.276 -21.260
## [92] -21.238 -21.214 -21.186 -21.158 -21.128 -21.099 -21.070
## [99] -21.043 -21.018
##
## $aar
## Time Series:
## Start = 4001
## End = 4100
## Frequency = 1
## [1] -129.938 -145.714 -152.940 -153.120 -148.085 -139.543 -128.812
## [8] -116.779 -104.019 -91.004 -78.215 -66.178 -55.329 -45.884
## [15] -37.844 -31.063 -25.352 -20.536 -16.475 -13.062 -10.219
## [22] -7.889 -6.028 -4.602 -3.577 -2.921 -2.599 -2.574
## [29] -2.806 -3.253 -3.876 -4.635 -5.494 -6.421 -7.388
## [36] -8.372 -9.355 -10.322 -11.262 -12.167 -13.034 -13.857
## [43] -14.638 -15.374 -16.068 -16.721 -17.335 -17.912 -18.455
## [50] -18.966 -19.448 -19.902 -20.332 -20.740 -21.126 -21.494
## [57] -21.845 -22.180 -22.501 -22.809 -23.105 -23.390 -23.666
## [64] -23.932 -24.189 -24.439 -24.681 -24.917 -25.147 -25.370
## [71] -25.588 -25.801 -26.009 -26.212 -26.411 -26.605 -26.795
## [78] -26.981 -27.163 -27.341 -27.516 -27.687 -27.854 -28.018
## [85] -28.179 -28.336 -28.490 -28.640 -28.788 -28.932 -29.072
## [92] -29.210 -29.345 -29.476 -29.604 -29.729 -29.852 -29.971
## [99] -30.087 -30.200
Figure 1. Forecasting Results for Different Models
plot(x.test)
lines(frc.test$linear, col = "red")
lines(frc.test$setar, col = "blue")
# lines(frc.test$lstar, col='green')
lines(frc.test$nnet, col = "orange")
lines(frc.test$aar, col = "violet")
Step 8. Construction of the Forecast Metrics.
sapply(mod.test, MAPE)
## linear setar nnet aar
## 0.3074 0.3029 0.3072 0.3069
rmse(x.test, frc.test$linear)
## [1] 93.45
rmse(x.test, frc.test$setar)
## [1] 86.36
# rmse(x.test,frc.test$lstar)
rmse(x.test, frc.test$nnet)
## [1] 94.55
rmse(x.test, frc.test$aar)
## [1] 94.09
Work in development at The Cromwell Workshop - www.cromwellworkshop.com by Jeff B. Cromwell, PhD