[R] StructTS with 2 seasons

Birgit Erni Birgit.Erni at uct.ac.za
Fri Dec 3 10:22:53 CET 2010


Thanks!

Following Prof. Ripley's reply I inserted the following lines into
StructTS (makeBSM)

if (nf < 3) 
           { ind <- 2L }
        else 
           { ind <- 3:nf
             T[cbind(ind + 1L, ind)] <- 1
           }

This worked and gave me the same variance estimates as the dlm package.


I also see that the StructTS function was updated and can handle 2
seasons now. 

Thanks for the help!

Birgit




I have replreplied that StructI have replaced t to the help from
Giovanni Petris and Brian Ripley


From: 	Giovanni Petris <gpetris at uark.edu>
To:	Birgit Erni <Birgit.Erni at uct.ac.za>
CC:	<r-help at r-project.org>
Date: 	30 November 2010 05:40 PM
Subject: 	Re: [R] StructTS with 2 seasons

Hi Birgit,

You are right, StructTS does not seem to be able to handle a time
series
with frequency 2 (and after a quick glance at the code I suspect a
frequency 3 may not work either).

There are contributed packages that you can use to fit (and smooth,
filter, forecast) structural time series and more general state space
models. Among them I would look at dlm, KFAS, dse, sspir.

HTH,
Giovanni
(author of 'dlm')

On Tue, 2010-11-30 at 10:40 +0200, Birgit Erni wrote:
> Dear All, 
> 
> I am trying to fit a structural time series model using the StructTS
> function (package stats) with only 2 seasons (summer and winter).
More
> than 2 seasons work fine but with 2 seasons I get this error:
> 
> > fit <- StructTS(y.ts, type="BSM")
> Error in T[cbind(ind + 1L, ind)] <- 1 : subscript out of bounds
> 
> I have looked at Prof. Ripley's 2002 RNews article but cannot see
why
> theoretically it should not be possible to have only 2 seasons. I
have
> also looked at the StructTS code, and the problem is with the
transition
> matrix in makeBSM, but I am not sure whether I should overwrite this
or
> whether my thinking is wrong. 
> 
> I attach code that simulates data for the basic structural model
(with
> time-varying trend and season), the top line s determines the number
of
> seasons in a year (3 and 4 work fine, 2 gives the above error), and
at
> the end tries to fit the model. I am using R version 2.11.1, Windows
XP.
> 
> 
> # ----- Basic structural model -----
> 
> s <- 2            # frequency or period
> 
> mu0 <- 6
> beta0 <- 1
> it0 <- 0
> it.states <- sample(1:6,s)
> 
> mu <- numeric()              # mean 
> beta <- numeric()           # slope
> it <- numeric()                # seasonal state
> 
> mu[1] <- mu0 + beta0 + rnorm(1,0,2)
> beta[1] <- beta0 + rnorm(1,0,0.5)
> it[1] <- - sum(it.states[2:s]) + rnorm(1,0,1)
> 
> it.states <- c(it[1],it.states[2:s])
> 
> for (i in 2:30)
>  { mu[i] <- mu[i-1] + beta[i-1] + rnorm(1,0,2)
>    beta[i] <- beta[i-1] + rnorm(1,0,2)
>    it[i] <- - sum(it.states[1:(s-1)]) + rnorm(1,0,1)
>    it.states <- c(it[i],it.states[1:(s-1)])
>  }
> 
> y <- numeric()
> 
> for (i in 1:30)
>  { y[i] <- mu[i] + it[i] + rnorm(1,0,3) }
> 
> par(mfrow=c(2,2))
> plot(y,type="b")
> plot(mu)
> plot(beta)
> plot(it)
> 
> y.ts <- ts(y,start=c(1977,1),frequency=s)
> y.ts
> 
> fit <- StructTS(y.ts, type="BSM")
> fit
> 
> # ------------------------------------------------------
> 
> Thank you.
> 
> Birgit Erni
> 
> Department of Statistical Sciences
> University of Cape Town, South Africa
> 
> 
> 
>  
> 
> ###
> UNIVERSITY OF CAPE TOWN 
> 
> This e-mail is subject to the UCT ICT policies and
e-mai...{{dropped:13}}
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help 
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html 
> and provide commented, minimal, self-contained, reproducible code.

-- 

Giovanni Petris  <GPetris at uark.edu>
Associate Professor
Department of Mathematical Sciences
University of Arkansas - Fayetteville, AR 72701
Ph: (479) 575-6324, 575-8630 (fax)
http://definetti.uark.edu/~gpetris/ 



 

###
UNIVERSITY OF CAPE TOWN 

This e-mail is subject to the UCT ICT policies and e-mai...{{dropped:13}}



More information about the R-help mailing list