R Notebook

A look at various filter and smoothing models in R with libraries from mass(structT, tsSmooth, kalmanLike etc.),by Venables & Ripley. KFAS, Kalman Filter and Smoother for Exponential Family State Space Models,by Jounie Helske,and atsa from SHumway and Stoffer to accompany their book Data sets and scripts to accompany Time Series Analysis and Its Applications: With R Examples (4th ed), by R.H. Shumway and D.S. Stoffer. Springer Texts in Statisics. Rob Hyndmans forecast package.(now updated as fable in tidymodels)

[1] "^VIX"

fitvix <- StructTS(vix, "level")
fitvix

Call:
StructTS(x = vix, type = "level")

Variances:
  level  epsilon  
  2.413    1.164  

StructTS(x = vix, type = "level")

Call:
StructTS(x = vix, type = "level")

Variances:
  level  epsilon  
  2.413    1.164  

plot(vix, col = "blue")
lines(fitted(fitvix), lty ="dashed", lwd =2)

lines(tsSmooth(fitvix),lty = "dotted", lwd =2)

plot(forecast(fitvix))

# rearranged not reason for kalman est unfit ie. forecast command

#lines(kalman_est, lty = "dotted", ,lwd = 2, col = "yellow")
# the kalman est fit did not run, syntax or logic error, for reader to pursue, kalman fit below.

# an ARIMA fit
fit3 <- arima(vix, c(1, 1, 0))
predict(fit3, 12)

$pred
Time Series:
Start = 2022.8796142289 
End = 2022.90973120598 
Frequency = 365.2425 
 [1] 24.02898 24.02134 24.02336 24.02283
 [5] 24.02297 24.02293 24.02294 24.02294
 [9] 24.02294 24.02294 24.02294 24.02294

$se
Time Series:
Start = 2022.8796142289 
End = 2022.90973120598 
Frequency = 365.2425 
 [1] 2.096842 2.604237 3.104520 3.516523
 [5] 3.889381 4.228451 4.542537 4.836204
 [9] 5.113048 5.375649 5.626008 5.865690

## reconstruct this
pr <- KalmanForecast(12, fit3$model)
pr$pred + fit3$coef[4]

 [1] NA NA NA NA NA NA NA NA NA NA NA NA

sqrt(pr$var * fit3$sigma2)

 [1] 2.096842 2.604237 3.104520 3.516523
 [5] 3.889381 4.228451 4.542537 4.836204
 [9] 5.113048 5.375649 5.626008 5.865690

## and now do it year by year
mod <- fit3$model
for(y in 1:3) {
  pr <- KalmanForecast(4, mod, TRUE)
  print(list(pred = pr$pred + fit3$coef["intercept"],
             se = sqrt(pr$var * fit3$sigma2)))
  mod <- attr(pr, "mod")
}

$pred
[1] NA NA NA NA

$se
[1] 2.096842 2.604237 3.104520 3.516523

$pred
[1] NA NA NA NA

$se
[1] 3.889381 4.228451 4.542537 4.836204

$pred
[1] NA NA NA NA

$se
[1] 5.113048 5.375649 5.626008 5.865690

fit_kalman <- KalmanLike(vix, mod, nit = 0L, update = FALSE)
out <- KalmanSmooth(vix, mod, nit = 0L)
kalman_est <- out$smooth[ ,2]
plot.ts(vix)

plot(kalman_est, lty = "dashed", lwd = 2,col = "green")

vix <-  ts_timeSeries(vix)
histPlot(vix)

library(KFAS)
library(tseries)
library(timeSeries)
library(zoo)
library(quantmod)
library(forecast)

Start Kalman code with structural time series. Using KFAS library

getSymbols("SPY",  from = as.Date("2021-01-01"),
           end = as.Date("2022-10-31"),
           quote = "Close")

[1] "SPY"

cls <- SPY$SPY.Close



sp_r <- dailyReturn(SPY)

sp_r <- as.ts(sp_r)

fit_structural  <-  SSModel(sp_r ~ SSMtrend(1, Q = list(matrix(NA))), H = matrix(NA))

 fit_structural <- fitSSM(fit_structural, inits = c(0, 0),method = "BFGS")

 fit_structural$model["Q"]

, , 1

             [,1]
[1,] 3.508162e-10

Plot returns model and start price series model and plot.

plot(fit_structural$model)

NA

fit_structural2  <-  SSModel(cls ~ SSMtrend(1, Q = list(matrix(NA))), H = matrix(NA))

Warning: setting this dimension may lead to an invalid zoo objectWarning: setting this dimension may lead to an invalid zoo object

fit_structural2 <- fitSSM(fit_structural2, inits = c(0, 0),method = "BFGS")

fit_structural2$model["Q"]

, , 1

             [,1]
[1,] 5.191219e-12

str(fit_structural2$model$y)

 Time-Series [1:474, 1] from 1 to 474: 369 371 374 379 381 ...
 - attr(*, "index")= num [1:474] 1.61e+09 1.61e+09 1.61e+09 1.61e+09 1.61e+09 ...
  ..- attr(*, "tzone")= chr "UTC"
  ..- attr(*, "tclass")= chr "Date"
 - attr(*, "src")= chr "yahoo"
 - attr(*, "updated")= POSIXct[1:1], format:  ...

plot(fit_structural2$model$y)

plot(fit_structural2$model$y, cls)

plot((fit_structural2$model$y), lty = "dotted", lwd =2)

fit <- forecast(fit_structural2$model$y, h = 10 )
plot(fit)

kalman filter using atsa library.

library(quantmod)
library(tidyverse)
library(tsbox)
library(tidyquant)
library(forecast)

start downloading RF data stock ticker Regional Finance corp. base comparison naive model forecast.

rf_price <- tq_get("RF", from = "2020-01-01")

rf_ts <- as.ts(rf_price)



n_fc <- naive(
  rf_price$close,
  h = 10,
  level = c(80, 95),
  fan = FALSE,
  lambda = NULL,
  biasadj = FALSE,
  x = rf_price$close
)

autoplot(n_fc)

NA
NA
NA

start structts code

  fit <- StructTS(rf_price$close, type = "level")
#if (struct$code != 0) stop("optimizer did not converge")

smoothed <- tsSmooth(fit)
str(smoothed)

 Time-Series [1:727, 1] from 1 to 727: 17.1 16.8 16.6 16.4 16.5 ...
 - attr(*, "dimnames")=List of 2
  ..$ : NULL
  ..$ : chr "level"

#tsdiag(fit)

plot(rf_price$close, type = "l" , lwd =1)
lines(fitted(fit), lty = "dashed", lwd = 2, col = "green")

lines(tsSmooth(fit), lty ="dotted", lwd = 2, col = "blue")

plot(forecast(fit, h= 20))

NA
NA

start atsa code

library(astsa)

fit2 <- ssm(rf_price$close,A=1, alpha=.01, phi=1, sigw=.01, sigv=.1)

initial  value 34062.188799 
iter   2 value 3019.274455
iter   3 value 2739.621359
iter   4 value 2734.440614
iter   5 value 2734.306521
iter   6 value 2734.285794
iter   7 value 2734.251764
iter   8 value 2734.159321
iter   9 value 2734.020029
iter  10 value 2733.891035
iter  11 value 2431.772512
iter  12 value 2364.036048
iter  13 value 2244.024934
iter  14 value 2037.962713
iter  15 value 1888.963511
iter  16 value 1763.906884
iter  17 value 1704.698280
iter  18 value 1590.163090
iter  19 value 1502.535715
iter  20 value 1385.516543
iter  21 value 1194.683236
iter  22 value 1178.227716
iter  23 value 905.085614
iter  24 value 505.799424
iter  25 value 272.220675
iter  26 value 148.561118
iter  27 value 49.117765
iter  28 value -14.156880
iter  29 value -49.523492
iter  30 value -57.536101
iter  31 value -125.439876
iter  32 value -163.507674
iter  33 value -171.330783
iter  34 value -171.696825
iter  35 value -175.095062
iter  36 value -175.249986
iter  37 value -177.864860
iter  38 value -179.788364
iter  39 value -179.799725
iter  40 value -179.805357
iter  41 value -179.805874
iter  42 value -179.809292
iter  43 value -180.007659
iter  44 value -180.011473
iter  45 value -180.053510
iter  46 value -180.075863
iter  47 value -180.077005
iter  48 value -180.077169
iter  48 value -180.077170
iter  48 value -180.077170
final  value -180.077170 
converged
        estimate          SE
phi   0.99555157 0.003689686
alpha 0.08899968 0.069553448
sigw  0.45969564 0.021371139
sigv  0.08071906 0.052033219

tsplot(rf_price$close, col=4, type="o", pch=20, ylab="RF Price")
lines(fit2$Xs, col=6, lwd=2)

xx = c(time(fit2$Xs), rev(time(fit2$Xs)))
yy = c(fit2$Xs-2*sqrt(fit2$Ps), rev(fit2$Xs+2*sqrt(fit2$Ps)))
polygon(xx, yy, border=8, col=gray(.6, alpha=.25) )

num = COUNT(rf_price$close)


Linn  = function(para){
  phi = para[1]; sigw = para[2]; sigv = para[3]
  kf  = Kfilter0(num, rf_price$close, A=1, mu0=0, Sigma0=10, phi, sigw, sigv)
  return(kf$like)
}


# Estimation
init.par = c(.1, 1, 1)   # initial parameter values
(est = optim(init.par, Linn, gr=NULL, method="BFGS", hessian=TRUE, control=list(trace=1,REPORT=1)))

initial  value 58153.511326 
iter   2 value 8985.289029
iter   3 value 8983.344225
iter   4 value 8978.601561
iter   5 value 8973.824843
iter   6 value 8969.013554
iter   7 value 8964.167153
iter   8 value 8959.285064
iter   9 value 8954.366750
iter  10 value 8949.411630
iter  11 value 8944.419064
iter  12 value 8939.388464
iter  13 value 8934.319194
iter  14 value 8929.210628
iter  15 value 8924.062088
iter  16 value 8918.872927
iter  17 value 8913.642421
iter  18 value 8908.369858
iter  19 value 8903.054484
iter  20 value 8897.695571
iter  21 value 8892.292344
iter  22 value 8886.843985
iter  23 value 8881.349670
iter  24 value 8875.808544
iter  25 value 8870.219782
iter  26 value 8864.582429
iter  27 value 8858.895609
iter  28 value 8853.158352
iter  29 value 8847.369672
iter  30 value 8841.528567
iter  31 value 8835.633991
iter  32 value 8829.684863
iter  33 value 8823.680077
iter  34 value 8817.618500
iter  35 value 8811.498926
iter  36 value 8805.320124
iter  37 value 8799.080851
iter  38 value 8792.779782
iter  39 value 8786.415595
iter  40 value 8779.986869
iter  41 value 8773.492137
iter  42 value 8766.929922
iter  43 value 8760.298683
iter  44 value 8753.596784
iter  45 value 8746.822564
iter  46 value 8739.974291
iter  47 value 8733.050150
iter  48 value 8726.048292
iter  49 value 8718.966728
iter  50 value 8711.803474
iter  51 value 8704.556422
iter  52 value 8697.223378
iter  53 value 8689.802046
iter  54 value 8682.290050
iter  55 value 8674.684908
iter  56 value 8666.984045
iter  57 value 8659.184738
iter  58 value 8651.284159
iter  59 value 8643.279363
iter  60 value 8635.167251
iter  61 value 8626.944575
iter  62 value 8618.607931
iter  63 value 8610.153750
iter  64 value 8601.578323
iter  65 value 8592.877671
iter  66 value 8584.047690
iter  67 value 8575.084010
iter  68 value 8565.982016
iter  69 value 8556.736942
iter  70 value 8547.343634
iter  71 value 8537.796673
iter  72 value 8528.090352
iter  73 value 8518.218621
iter  74 value 8508.175046
iter  75 value 8497.952836
iter  76 value 8487.544709
iter  77 value 8476.942938
iter  78 value 8466.139298
iter  79 value 8455.125018
iter  80 value 8443.890691
iter  81 value 8432.426232
iter  82 value 8420.720826
iter  83 value 8408.762901
iter  84 value 8396.539923
iter  85 value 8384.038436
iter  86 value 8371.243830
iter  87 value 8358.140323
iter  88 value 8344.710776
iter  89 value 8330.936533
iter  90 value 8316.797196
iter  91 value 8302.270522
iter  92 value 8287.332082
iter  93 value 8271.954993
iter  94 value 8256.109582
iter  95 value 8239.762962
iter  96 value 8222.878669
iter  97 value 8205.415901
iter  98 value 8187.328996
iter  99 value 8168.566582
iter 100 value 8149.070508
final  value 8149.070508 
stopped after 100 iterations
$par
[1] 426.1884 510.7464 173.0802

$value
[1] 8149.071

$counts
function gradient 
     103      100 

$convergence
[1] 1

$message
NULL

$hessian
              [,1]         [,2]
[1,] -1.232593e-03 1.397552e-03
[2,]  1.397552e-03 2.807610e-03
[3,] -2.262368e-05 8.003553e-05
              [,3]
[1,] -2.262368e-05
[2,]  8.003553e-05
[3,] -2.170464e-02

SE = sqrt(diag(solve(est$hessian)))

Warning: NaNs produced

cbind(estimate=c(phi=est$par[1],sigw=est$par[2],sigv=est$par[3]), SE)

     estimate       SE
phi  426.1884      NaN
sigw 510.7464 15.08872
sigv 173.0802      NaN

ks =  Ksmooth0(num, rf_price$close, A=1,mu0=0, Sigma0=10, est$par[1], est$par[2], est$par[3])
tsplot(rf_price$close, type='o', col=4, pch=20)
lines(ks$xf, col=6, lwd=2)

NA
NA
NA
NA
NA
NA