helske / KFAS

KFAS: R Package for Exponential Family State Space Models
64 stars 17 forks source link

SSMSeasonal harmonics wrong indices #30

Closed jyoder6 closed 6 years ago

jyoder6 commented 6 years ago

I believe the harmonics indexing logic is bugged:


library(KFAS)
seas.full <- SSMseasonal(12, Q=.1, sea.type='trigonometric', P1=1e6)
seas.full$state_names
# [1] "sea_trig1"  "sea_trig*1" "sea_trig2"  "sea_trig*2" "sea_trig3"  "sea_trig*3" "sea_trig4"  "sea_trig*4"
# [9] "sea_trig5"  "sea_trig*5" "sea_trig6" 
# Suppose I want the 1st, 2nd, and 4th pairs of functions to get the slowest, next slowest, and 4th slowest frequencies.
harmonics = c(1, 2, 4)
# Wanted indices == 
names.wanted = c("sea_trig1",  "sea_trig*1", "sea_trig2",  "sea_trig*2", "sea_trig4",  "sea_trig*4")
ind.wanted = which(seas.full$state_names %in% names.wanted)
# [1] 1 2 3 4 7 8

seas.actual <- SSMseasonal(12, harmonics=harmonics, Q=.1, sea.type='trigonometric', P1=1e6)
# Note this warning:
# Warning message:
  # In rep(harmonics, each = 2) * rep(1:2, each = 2) :
  # longer object length is not a multiple of shorter object length
ind.actual <- which(seas.actual$state_names %in% names.wanted)

assertthat::are_equal(ind.actual, ind.wanted)
# [1] FALSE

# Double check:
seas.actual$state_names
# [1] "sea_trig1"  "sea_trig*1" "sea_trig*2" "sea_trig3"  "sea_trig*2" "sea_trig3" 

# Looking at the source code
# in SSMseasonal.R line 142
# is
ind.actual.142 <- rep(harmonics, each = 2) * rep(1:2, each = 2) + 0:1
# Note warning message here

# Should be
ind.shouldbe.142 <- rep(harmonics, each = 2) * 2 - 1:0

assertthat::are_equal(ind.actual.142, ind.wanted)
# [1] FALSE
assertthat::are_equal(ind.shouldbe.142, ind.wanted)
# [1] TRUE

This relates to #21. I am OK with making a pull request but didn't see contributing instructions so I am putting this in as an issue.

helske commented 6 years ago

Thanks, the original implementation was quite a quick solution without much testing, that rep line is obviously wrong. Feel free to make a pull request, thank you in advance!