The SSM Toolbox for Matlab

The SSM T oolbox for Matlab Jyh-Y ing Peng 1 , 2 , John A. D. Aston 1 1 Institute of Statistical Science, Academia Sinica ∗ 2 Computer Science and Information Engineering, National T aiwan Uni versity Nov ember 26, 2024 ∗ Address for Correspondence: Institute of Statistical Science, Academia Sinica 128 Academia Road, Sec 2 T aipei 115, T aiwan email: { jypeng,jaston } @stat.sinica.edu.tw 2 Peng & Aston Contents 1 Introduction 5 2 The State Space Models 5 2.1 Linear Gaussian models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2 Non-Gaussian models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.3 Nonlinear models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 3 The State Space Algorithms 6 3.1 Kalman ﬁlter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 3.2 State smoother . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 3.3 Disturbance smoother . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 3.4 Simulation smoother . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 4 Getting Started 10 4.1 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 4.2 Model construction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 4.3 Model and state estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 5 SSM Classes 11 5.1 Class SSMA T . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 5.2 Class SSDIST . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 5.3 Class SSFUNC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 5.4 Class SSP ARAM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 5.5 Class SSMODEL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 6 SSM Functions 14 7 Data Analysis Examples 15 7.1 Structural time series models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 7.1.1 Univ ariate analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 7.1.2 Biv ariate analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 7.2 ARMA models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 7.3 Cubic spline smoothing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 7.4 Poisson distribution error model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 7.5 t-distribution models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 7.6 Hillmer-T iao decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 A Predeﬁned Model Reference 29 B Class Reference 32 B.1 SSMA T . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 B.1.1 Subscripted reference and assignment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 B.1.2 Class methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 B.2 SSDIST . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 B.2.1 Subscripted reference and assignment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 SSM T oolbox for Matlab 3 B.2.2 Class methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 B.2.3 Predeﬁned non-Gaussian distrib utions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 B.3 SSFUNC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 B.3.1 Subscripted reference and assignment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 B.3.2 Class methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 B.4 SSP ARAM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 B.4.1 Subscripted reference and assignment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 B.4.2 Class methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 B.5 SSMODEL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 B.5.1 Subscripted reference and assignment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 B.5.2 Class Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 C Function Reference 38 C.1 User-deﬁned update functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 C.2 Data analysis functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 C.3 Stock element functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 C.4 Miscellaneous helper functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4 Peng & Aston Abstract State Space Models (SSM) is a MA TLAB 7.0 software toolbox for doing time series analysis by state space meth- ods. The software features fully interactiv e construction and combination of models, with support for univariate and multiv ariate models, complex time-varying (dynamic) models, non-Gaussian models, and various standard models such as ARIMA and structural time-series models. The software includes standard functions for Kalman ﬁltering and smoothing, simulation smoothing, lik elihood ev aluation, parameter estimation, signal extraction and forecasting, with incorporation of exact initialization for ﬁlters and smoothers, and support for missing observations and multiple time series input with common analysis structure. The software also includes implementations of TRAMO model selection and Hillmer-T iao decomposition for ARIMA models. The softw are will provide a general toolbox for doing time series analysis on the MA TLAB platform, allo wing users to tak e adv antage of its readily av ailable graph plotting and general matrix computation capabilities. Keyw ords : T ime Series Analysis, State Space Models, Kalman Filter . SSM T oolbox for Matlab 5 1 Intr oduction State Space Models (SSM) is a MA TLAB toolbox for doing time series analysis using general state space models and the Kalman ﬁlter [1]. The goal of this software package is to provide users with an intuitive, conv enient and efﬁcient way to do general time series modeling within the state space framework. Speciﬁcally , it seeks to provide users with easy construction and combination of arbitrary models without having to e xplicitly deﬁne ev ery component of the model, and to maximize transparency in their data analysis usage so no special consideration is needed for any individual model. This is achiev ed through the uniﬁcation of all state space models and their extension to non- Gaussian and nonlinear special cases. The user creation of custom models is also implemented to be as general, ﬂexible and efﬁcient as possible. Thus, there are often multiple ways of deﬁning a single model and choices as to the parametrization versus initialization and to how the model update functions are implemented. Stock model components are also provided to ease user extension to e xisting predeﬁned models. Functions that implement standard algorithms such as the Kalman ﬁlter and state smoother , log likelihood calculation and parameter estimation will work seamlessly across all models, including any user deﬁned custom models. These features are implemented through object-oriented programming primiti ves provided by MA TLAB and classes in the toolbox are deﬁned to conform to MA TLAB conv entions whenev er possible. The result is a highly integrated toolbox with support for general state space models and standard state space algorithms, complemented by the built-in matrix computation and graphic plotting capabilities of MA TLAB. 2 The State Space Models This section presents a summary of the basic deﬁnition of models supported by SSM. Currently SSM implements the Kalman ﬁlter and related algorithms for model and state estimation, hence non-Gaussian or nonlinear models need to be approximated by linear Gaussian models prior to or during estimation. Howe ver the approximation is done automatically and seamlessly by the respectiv e routines, e ven for user -deﬁned non-Gaussian or nonlinear models. The following notation for v arious sequences will be used throughout the paper: y t p × 1 Observation sequence ε t p × 1 Observation disturbance (Unobserv ed) α t m × 1 State sequence (Unobserv ed) η t r × 1 State disturbance (Unobserved) 2.1 Linear Gaussian models SSM supports linear Gaussian state space model in the form y t = Z t α t + ε t , ε t ∼ N (0 , H t ) , α t +1 = c t + T t α t + R t η t , η t ∼ N (0 , Q t ) , α 1 ∼ N ( a 1 , P 1 ) , t = 1 , . . . , n. (1) Thus the matrices Z t , c t , T t , R t , H t , Q t , a 1 , P 1 are required to deﬁne a linear Gaussian state space model. The matrix Z t is the state to observation linear transformation, for univ ariate models it is a ro w vector . The matrix c t is the same size as the state vector , and is the “constant” in the state update equation, although it can be dynamic or dependent on model parameters. The square matrix T t deﬁnes the time ev olution of states. The matrix R t transforms general disturbance into state space, and exists to allow for more varieties of models. H t and Q t are Gaussian variance matrices gov erning the disturbances, and a 1 and P 1 are the initial conditions. The matrices and their dimensions are listed: 6 Peng & Aston Z t p × m State to observation transform matrix c t m × 1 State update constant T t m × m State update transform matrix R t m × r State disturbance transform matrix H t p × p Observation disturbance v ariance Q t r × r State disturbance variance a 1 m × 1 Initial state mean P 1 m × m Initial state v ariance 2.2 Non-Gaussian models SSM supports non-Gaussian state space models in the form y t ∼ p ( y t | θ t ) , θ t = Z t α t , α t +1 = c t + T t α t + R t η t , η t ∼ Q t = p ( η t ) , α 1 ∼ N ( a 1 , P 1 ) , t = 1 , . . . , n. (2) The sequence θ t is the signal and Q t is a non-Gaussian distribution (e.g. heavy-tailed distribution). The non-Gaussian observation disturbance can tak e two forms: an exponential family distrib ution H t = p ( y t | θ t ) = exp  y T t θ t − b t ( θ t ) + c t ( y t )  , −∞ < θ t < ∞ , (3) or a non-Gaussian additiv e noise y t = θ t + ε t , ε t ∼ H t = p ( ε t ) . (4) W ith model combination it is also possible for H t and Q t to be a combination of Gaussian distributions (represented by variance matrices) and v arious non-Gaussian distributions. 2.3 Nonlinear models SSM supports nonlinear state space models in the form y t = Z t ( α t ) + ε t , ε t ∼ N (0 , H t ) , α t +1 = c t + T t ( α t ) + R t η t , η t ∼ N (0 , Q t ) , α 1 ∼ N ( a 1 , P 1 ) , t = 1 , . . . , n. (5) Z t and T t are functions that map m × 1 vectors to p × 1 and m × 1 vectors respectiv ely . W ith model combination it is also possible for Z t and T t to be a combination of linear functions (matrices) and nonlinear functions. 3 The State Space Algorithms This section summarizes the core algorithms implemented in SSM, starting with the Kalman ﬁlter . Normally the initial conditions a 1 and P 1 should be given for these algorithms, but the exact diffuse initialization implemented permits elements of either to be unknown, and deri ved implicitly from the ﬁrst few observations. The unkno wn elements in a 1 are set to 0 , and the corresponding diagonal entry in P 1 is set to ∞ , this represents an improper prior which reﬂects the lack of kno wledge about that particular element a priori . See [1] for details. T o keep v ariables ﬁnite in the algorithms we deﬁne another initial condition P ∞ , 1 with elements ( P ∞ , 1 ) i,j = ( 0 , if ( P 1 ) i,j < ∞ , 1 , if ( P 1 ) i,j = ∞ , (6) SSM T oolbox for Matlab 7 and set inﬁnite elements of P 1 to 0 . All the state space algorithms supports exact diffuse initialization. For non-Gaussian and nonlinear models linear Gaussian approximations are performed by linearization of the loglikelihood equation and model matrices respectiv ely , and the Kalman ﬁlter algorithms can then be used on the resulting approximation models. 3.1 Kalman ﬁlter Giv en the observation y t , the model and initial conditions a 1 , P 1 , P ∞ , 1 , Kalman ﬁlter outputs a t = E( α t | y 1: t − 1 ) and P t = V ar( α t | y 1: t − 1 ) , the one step ahead state prediction and prediction variance. Theoretically the dif fuse state variance P ∞ ,t will only be non-zero for the ﬁrst couple of time points (usually the number of diffuse elements, b ut may be longer depending on the structure of Z t ), after which exact diffuse initialization is completed and normal Kalman ﬁlter is applied. Here n denotes the length of observation data, and d denotes the number of time points where the state variance remains dif fuse. The normal Kalman ﬁlter recursion (for t > d ) are as follows: v t = y t − Z t a t , F t = Z t P t Z T t + H t , K t = T t P t Z T t F − 1 t , L t = T t − K t Z t , a t +1 = c t + T t a t + K t v t , P t +1 = T t P t L T t + R t Q t R T t . The exact dif fuse initialized Kalman ﬁlter recursion (for t ≤ d ) is divided into two cases, when F ∞ ,t = Z t P ∞ ,t Z T t is nonsingular: v t = y t − Z t a t , F (1) t = ( Z t P ∞ ,t Z T t ) − 1 , F (2) t = − F (1) t ( Z t P t Z T t + H t ) F (1) t , K t = T t P ∞ ,t Z T t F (1) t , K (1) t = T t ( P t Z T t F (1) t + P ∞ ,t Z T t F (2) t ) , L t = T t − K t Z t , L (1) t = − K (1) t Z t , a t +1 = c t + T t a t + K t v t , P t +1 = T t P ∞ ,t L (1) T t + T t P t L T t + R t Q t R T t , P ∞ ,t +1 = T t P ∞ ,t L T t , and when F ∞ ,t = 0 : v t = y t − Z t a t , F t = Z t P t Z T t + H t , K t = T t P t Z T t F − 1 t , L t = T t − K t Z t , a t +1 = c t + T t a t + K t v t , 8 Peng & Aston P t +1 = T t P t L T t + R t Q t R T t , P ∞ ,t +1 = T t P ∞ ,t T T t . 3.2 State smoother Giv en the Kalman ﬁlter outputs, the state smoother outputs ˆ α t = E( α t | y 1: n ) and V t = V ar( α t | y 1: n ) , the smoothed state mean and variance gi ven the whole observation data, by backw ards recursion. T o start the backwards recursion, ﬁrst set r n = 0 and N n = 0 , where r t and N t are the same size as a t and P t respectiv ely . The normal backwards recursion (for t > d ) is then r t − 1 = Z T t F − 1 t v t + L T t r t , N t − 1 = Z T t F − 1 t Z t + L T t N t L t , ˆ α t = a t + P t r t − 1 , V t = P t − P t N t − 1 P t . For the exact diffuse initialized portion ( t ≤ d ) set additional diffuse variables r (1) d = 0 and N (1) d = N (2) d = 0 corresponding to r d and N d respectiv ely , the backwards recursion is r (1) t − 1 = Z T t F (1) t v t + L T t r (1) t + L (1) T t r t , r t − 1 = L T t r t , N (2) t − 1 = Z T t F (2) t Z t + L T t N (2) t L t + L T t N (1) t L (1) t + L (1) T t N (1) t L t + L (1) T t N t L (1) t , N (1) t − 1 = Z T t F (1) t Z t + L T t N (1) t L t + L (1) T t N t L t , N t − 1 = L T t N t L t , ˆ α t = a t + P t r t − 1 + P ∞ ,t r (1) t − 1 , V t = P t − P t N t − 1 P t − ( P ∞ ,t N (1) t − 1 P t ) T − P ∞ ,t N (1) t − 1 P t − P ∞ ,t N (2) t − 1 P ∞ ,t , if F ∞ ,t is nonsingular , and r (1) t − 1 = T T t r (1) t , r t − 1 = Z T t F − 1 t v t + L T t r t , N (2) t − 1 = T T t N (2) t T t , N (1) t − 1 = T T t N (1) t L t , N t − 1 = Z T t F − 1 t Z t + L T t N t L t , ˆ α t = a t + P t r t − 1 + P ∞ ,t r (1) t − 1 , V t = P t − P t N t − 1 P t − ( P ∞ ,t N (1) t − 1 P t ) T − P ∞ ,t N (1) t − 1 P t − P ∞ ,t N (2) t − 1 P ∞ ,t , if F ∞ ,t = 0 . 3.3 Disturbance smoother Giv en the Kalman ﬁlter outputs, the disturbance smoother outputs ˆ ε t = E( ε t | y 1: n ) , V ar( ε t | y 1: n ) , ˆ η t = E( η t | y 1: n ) and V ar( η t | y 1: n ) . Calculate the sequence r t and N t by backwards recursion as before, the disturbance can then be SSM T oolbox for Matlab 9 estimated by ˆ η t = Q t R T t r t , V ar( η t | y 1: n ) = Q t − Q t R T t N t R t Q t , ˆ ε t = H t ( F − 1 t v t − K T t r t ) , V ar( ε t | y 1: n ) = H t − H t ( F − 1 t + K T t N t K t ) H t , if t > d or F ∞ ,t = 0 , and ˆ η t = Q t R T t r t , V ar( η t | y 1: n ) = Q t − Q t R T t N t R t Q t , ˆ ε t = − H t K T t r t , V ar( ε t | y 1: n ) = H t − H t K T t N t K t H t , if t ≤ d and F ∞ ,t is nonsingular . 3.4 Simulation smoother The simulation smoother randomly generates an observation sequence ˜ y t and the underlying state and disturbance sequences ˜ α t , ˜  t and ˜ η t conditional on the model and observ ation data y t . The algorithm is as follows (each step is performed for all time points t ): 1. Draw random samples of the state and observ ation disturbances: ε + t ∼ N (0 , H t ) , η + t ∼ N (0 , Q t ) . 2. Draw random samples of the initial state: α + 1 ∼ N ( a 1 , P 1 ) . 3. Generate state sequences α + t and observation sequences y + t from the random samples. 4. Use disturbance smoothing to obtain ˆ ε t = E( ε t | y 1: n ) , ˆ η t = E( η t | y 1: n ) , ˆ ε + t = E( ε + t | y + 1: n ) , ˆ η + t = E( η + t | y + 1: n ) . 5. Use state smoothing to obtain ˆ α t = E( α t | y 1: n ) , ˆ α + t = E( α + t | y + 1: n ) . 6. Calculate ˜ α t = ˆ α t − ˆ α + t + α + t , ˜ ε t = ˆ ε t − ˆ ε + t + ε + t , ˜ η t = ˆ η t − ˆ η + t + η + t . 7. Generate ˜ y t from the output. 10 Peng & Aston 4 Getting Started The easiest and most frequent way to start using SSM is by constructing predeﬁned models, as opposed to creating a model from scratch. This section presents some examples of simple time series analysis using predeﬁned models, the complete list of av ailable predeﬁned models can be found in appendix A. 4.1 Preliminaries Some parts of SSM is programmed in c for maximum efﬁciency , before using SSM, go to the subfolder src and ex ecute the script mexall to compile and distribute all needed c me x functions. 4.2 Model construction T o create an instance of a predeﬁned model, call the SSMODEL (state space model) constructor with a coded string representing the model as the ﬁrst argument, and optional ar guments as necessary . For e xample: • model = ssmodel(’llm’) creates a local level model. • model = ssmodel(’arma’, p, q) creates an ARMA ( p, q ) model. The resulting variable model is a SSMODEL object and can be displayed just like any other MA TLAB variables. T o set or change the model parameters, use model.param , which is a row vector that behav es like a MA TLAB matrix except its size cannot be changed. The initial conditions usually defaults to exact diffuse initialization, where model.a1 is zero, and model.P1 is ∞ on the diagonals, but can likewise be changed. Models can be combined by horizontal concatenation, where only the observ ation disturbance model of the ﬁrst one will be retained. The class SSMODEL will be presented in detail in later sections. 4.3 Model and state estimation W ith the model created, estimation can be performed. SSM expects the data y to be a matrix of dimensions p × n , where n is the data size (or time duration). The model parameters are estimated by maximum likelihood, the SSMODEL class method estimate performs the estimation. For e xample: • model1 = estimate(y, model0) estimates the model and stores the result in model1 , where the pa- rameter values of model0 is used as initial value. • [model1 logL] = estimate(y, model0, psi0, [], optname1, optvalue1, optname2, optvalue2, ...) estimates the model with psi0 as the initial parameters using option settings speciﬁed with option value pairs, and returns the resulting model and loglik elihood. After the model is estimated, state estimation can be performed, this is done by the SSMODEL class method kalman and statesmo , which is the Kalman ﬁlter and state smoother respectiv ely . • [a P] = kalman(y, model) applies the Kalman ﬁlter on y and returns the one-step-ahead state predic- tion and variance. • [alphahat V] = statesmo(y, model) applies the state smoother on y and returns the expected state mean and variance. The ﬁltered and smoothed state sequences a and alphahat are m × n+1 and m × n matrices respectiv ely , and the ﬁltered and smoothed state v ariance sequences P and V are m × m × n+1 and m × m × n matrices respectively , except if m = 1 , in which case they are squeezed and transposed. SSM T oolbox for Matlab 11 5 SSM Classes SSM is composed of v arious MA TLAB classes each of which represents speciﬁc parts of a general state space model. The class SSMODEL represents a state space model, and is the most important class of SSM. It embeds the classes SSMAT , SSDIST , SSFUNC and SSPARAM , all designed to facilitate the construction and usage of SSMODEL for general data analysis. Consult appendix B for complete details about these classes. 5.1 Class SSMA T The class state space matrix ( SSMAT ) forms the basic components of a state space model, they can represent linear transformations that map vectors to v ectors (the Z t , T t and R t matrices), or variance matrices of Gaussian disturbances (the H t and Q t matrices). The ﬁrst data member of class SSMAT is mat , which is a matrix that represents the stationary part of the state space matrix, and is often the only data member that is needed to represent a state space matrix. F or dynamic state space matrices, the additional data members dmmask and dvec are needed. dmmask is a logical matrix the same size as mat , and speciﬁes elements of the state space matrix that varies with time, these elements arranged in a column in column-wise order forms the columns of the matrix dvec , thus nnz(dmmask) is equal to the number of rows in dvec , and the number of columns of dvec is the number of time points speciﬁed for the state space matrix. This is designed so that the code mat(dmmask) = dvec(:, t); gets the value of the state space matrix at time point t . Now these are the constant part of the state space matrix, in that they do not depend on model parameters. The logical matrix mmask speciﬁes which elements of the stationary mat is dependent on model parameters (variable), while the logical column vector dvmask speciﬁes the variable elements of dvec . Functions that update the state space matrix take the model parameters as input argument, and should output results so that the following code correctly updates the state space matrix: mat(mmask) = func1(psi); dvec(dvmask, :) = func2(psi); where func1 update the stationary part, and func2 updates the dynamic part. The class SSMAT is designed to represent the “structure” and current value of a state space matrix, thus the update functions are not part of the class and will be described in detail in later sections. Objects of class SSMAT behaves like a matrix to a certain extent, calling size returns the size of the stationary part, concatenations horzcat , vertcat and blkdiag can be used to combine SSMAT with each other and ordi- nary matrices, and same for plus . Data members mat and dvec can be read and modiﬁed provided the structure of the state space matrix stays the same, parenthesis indexing can also be used. 5.2 Class SSDIST The class state space distributions ( SSDIST ) is a child class of SSMAT , and represents non-Gaussian distributions for the observ ation and state disturbances ( H t and Q t ). Since disturbances are v ectors, it’ s possible to have both elements with Gaussian distributions, and elements with different non-Gaussian distributions. Also in order to use the Kalman ﬁlter , non-Gaussian distributions are approximated by Gaussian noise with dynamic v ariance. Thus the class SSDIST needs to keep track of multiple non-Gaussian distrib utions and also their corresponding disturbance elements, so that the individual Gaussian approximations can be concatenated in the correct order . 12 Peng & Aston As non-Gaussian distrib utions are approximated by dynamic Gaussian distrib ution in SSM, its representation need only consist of data required in the approximation process, namely the functions that take disturbance samples as input and outputs the dynamic Gaussian variance. Also a function that outputs the probability of the disturbance samples is needed for loglikelihood calculation. These functions are stored in two data members matf and logpf , both cell arrays of function handles. A logical matrix diagmask is needed to keep track of which elements each non-Gaussian distribution gov erns, where each column of diagmask corresponds to each distribution. Elements not speciﬁed in any columns are therefore Gaussian. Lastly a logical vector dmask of the same length as matf and logpf speciﬁes which distributions are v ariable, for example, t-distrib utions can be variable if the variance and de gree of freedom are model parameters. The update function func for SSDIST objects should be deﬁned to take parameter values as input arguments and output a cell matrix of function handles such that distf = func(psi); [matf { dmask } ] = distf { :, 1 } ; [logpf { dmask } ] = distf { :, 2 } ; correctly updates the non-Gaussian distributions. In SSM most non-Gaussian distributions such as exponential family distributions and some heavy-tailed noise distributions are predeﬁned and can be constructed similarly to predeﬁned models, these can then be inserted in stock components to allow for non-Gaussian disturbance generalizations. 5.3 Class SSFUNC The class state space functions ( SSFUNC ) is another child class of SSMAT , and represents nonlinear functions that transform vectors to vectors ( Z t , T t and R t ). As with non-Gaussian distribution it is possible to “concatenate” non- linear functions corresponding to dif ferent elements of the state v ector and destination vector . The dimension of input and output vector and which elements they corresponds to will hav e to be kept track of for each nonlinear function, and their linear approximations are concatenated using these information. Both the function itself and its deriv ati ve are needed to calculate the linear approximation, these are stored in data members f and df , both cell arrays of function handles. horzmask and vertmask are both logical matrices where each column speciﬁes the elements of the input and output vector for each function, respectiv ely . In this way the dimensions of the approximating matrix for the i th function can be determined as nnz(vertmask(:, i)) × nnz(horzmask(:, i)) . The logical vector fmask speciﬁes which functions are variable with respect to model parameters. The update function func for SSFUNC objects should be deﬁned to take parameter values as input arguments and output a cell matrix of function handles such that funcf = func(psi); [f { fmask } ] = funcf { :, 1 } ; [df { fmask } ] = funcf { :, 2 } ; correctly updates the nonlinear functions. 5.4 Class SSP ARAM The class state space parameters ( SSPARAM ) stores and organizes the model parameters, each model parameter has a name, v alue, transform and possible constraints. Transforms are necessary because most model parameters hav e a restricted domain, but most numerical optimization routines generally operate over the whole real line. For example a variance parameter σ 2 are suppose to be positi ve, so if σ 2 is optimized directly there’ s the possibility of error, but SSM T oolbox for Matlab 13 by deﬁning ψ = log( σ 2 ) and optimizing ψ instead, this problem is avoided. Model parameter transforms are not restricted to single parameters, sometimes a group of parameters need to be transformed together , the same goes for constraints. Hence the class SSPARAM also keeps track of parameter grouping. Like SSMAT , SSPARAM objects can be treated like a ro w vector of parameters to a certain extent. Speciﬁcally , it is possible to horizontally concatenate SSPARAM objects as a means of combining parameter groups. Other operations on SSPARAM objects include getting and setting the parameter values and transformed values. The untransformed values are the original parameter v alues meaningful to the model deﬁned, and will be referred to as param in code examples; the transformed v alues are meant only as input to optimization routines and hav e no direct interpretation in terms of the model, they are referred to as psi in code examples, but the user will rarely need to use them directly if predeﬁned models are used. 5.5 Class SSMODEL The class state space model ( SSMODEL ) represents state space models by embedding SSMAT , SSDIST , SSFUNC and SSPARAM objects and also keeping track of update functions, parameter masks and general model component information. The basic constituent elements of a model is described earlier as Z , T , R , H , Q , c , a1 and P1 . Z , T and R are vector transformations, the ﬁrst two can be SSFUNC or SSMAT objects, but the last one can only be a SSMAT object. H and Q govern the disturbances and can be SSDIST or SSMAT objects. c , a1 and P1 can only be SSMAT objects. These form the “constant” part of the model speciﬁcation, and is the only part needed in Kalman ﬁlter and related state estimation algorithms. The second part of SSMODEL concerns the update functions and model parameters, general speciﬁcations for the model estimation process. func and grad are cell arrays of update functions and their gradient for the basic model elements respectively . It is possible for a single function to update multiple parts of the model and these information are stored in the data member A , a length(func) × 18 adjacency matrix. Each row of A represents one update function, and each column represents one updatable element. The 18 updatable model elements in order are H Z T R Q c a1 P1 Hd Zd Td Rd Qd cd Hng Qng Znl Tnl where the d suf ﬁx means dynamic, ng means non-Gaussian and nl means nonlinear . It is therefore possible for a function that updates Z to output [vec dsubvec funcf] updating the stationary part of Z , dynamic part of Z and the nonlinear functions in Z respectiv ely . The row in A for this function would be H Z T R Q c a1 P1 Hd Zd Td Rd Qd cd Hng Qng Znl Tnl 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 Any function can update any combination of these 18 elements, the only restriction is that the output should follow the order set in A with possible omissions. On the other hand A also allo ws many functions to update a single model element, this is the main mechanism that enables smooth model combination. Suppose T1 of model 1 is updated by func1 , and T2 of model 2 by func2 , combination of models 1 and 2 requires block diagonal concatenation of T1 and T2 : T = blkdiag(T1, T2) . It is easy to see that the vertical concatenation of the output of func1 and func2 correctly updates the new matrix T as follows vec1 = func1(psi); vec2 = func2(psi); T.mat(T.mmask) = [vec1; vec2]; 14 Peng & Aston This also holds for horizontal concatenation, but for vertical concatenation of more than one column this would fail, luckily this case ne ver occurs when combining state space models, so can be safely ignored. The updates of dynamic, non-Gaussian and nonlinear elements can also be combined in this way , in fact, the whole scheme of the update functions are designed with combinable output in mind. For an example adjacency matrix A H Z T R Q c a1 P1 Hd Zd Td Rd Qd cd Hng Qng Znl Tnl func1 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 func2 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 func3 0 1 1 0 1 0 0 0 0 1 0 0 0 0 0 0 1 0 model update would then proceed by ﬁrst in voking the three update functions [Hvec1 Tvec1 Qvec1 Zfuncf1] = func1(psi); [Tvec2 Zdsubvec2] = func2(psi); [Zvec3 Tvec3 Qvec3 Zdsubvec3 Zfuncf3] = func3(psi); and then update the various elements H.mat(H.mmask) = Hvec1; Z.mat(Z.mmask) = Zvec3; T.mat(T.mmask) = [Tvec1; Tvec2; Tvec3]; Q.mat(Q.mmask) = [Qvec1; Qvec3]; Z.dvec(Z.dvmask, :) = [Zdsubvec2; Zdsubvec3]; Z.f(Z.fmask) = [Zfuncf1 { :, 1 } ; Zfuncf3 { :, 1 } ]; Z.df(Z.fmask) = [Zfuncf1 { :, 2 } ; Zfuncf3 { :, 2 } ]; Now to facilitate model construction and alteration, functions that need adjacency matrix A as input actually expects another form based on strings. Each of the 18 elements shown earlier is represented by the corresponding strings, ’H’ , ’Qng’ , ’Zd’ , . . . etc., and each row of A is a single string in a cell array of strings, where each string is a concatenation of v arious elements in any order . This form (a cell array of strings) is called “adjacency string” and is used in all functions that require an adjacency matrix as input. The remaining data members of SSMODEL are psi and pmask . psi is a SSPARAM object that stores the model parameters, and pmask is a cell array of logical row v ectors that speciﬁes which parameters are required by each up- date function. So each update function are in fact provided a subset of model parameters func i (psi.value(pmask { i } )) . It is trivial to combine these tw o data members when combining models. In summary the class SSMODEL can be divided into two parts, a “constant” part that corresponds to an instance of a theoretical state space model giv en all model parameters, which is the only part required to run Kalman ﬁlter and related state estimation routines, and a “variable” part that keeps track of the model parameters and how each parameter effects the values of various model elements, which is used (in addition to the “constant” part) in model estimation routines that produce parameter estimates. This division is the careful separation of the “structure” of a model and its “instantiation” giv en model parameters, and allows logical separation of various structural elements of the model (with the classes SSMAT , SSDIST and SSFUNC ) to facilitate model combination and usage, without compromising the ability to efﬁciently deﬁne and combine comple x model update mechanisms. 6 SSM Functions There are several types of functions in SSM besides class methods, these are data analysis functions, stock element functions and various helper functions. Most data analysis functions are Kalman ﬁlter related and implemented as SSM T oolbox for Matlab 15 SSMODEL methods, but for the user this is transparent since class methods and functions are called in the same way in MA TLAB. Examples for data analysis functions include kalman , statesmo and loglik , all of which uses the kalman ﬁlter and outputs the ﬁltered states, smoothed states and loglikelihood respectively . These and similar functions take three arguments, the data y , which must be p × n , where p is the number of variables in the data and n the data size, the model model , a SSMODEL object, and optional analysis option name value pairs (see function setopt in Appendix C.4), where settings such as tolerance can be speciﬁed. The functions linear and gauss calculates the linear Gaussian approximation models, and requires y , model and optional initial state estimate alpha0 as input arguments. The function estimate runs numerical optimization routines for model estimation; model updates, necessary approximations and loglikelihood calculations are done automatically in the procedure, input arguments are y , model and initial parameter estimate psi0 , the function outputs the maximum likelihood model. There are also some functions speciﬁc to ARIMA-type models, such as Hillmer-T iao decomposition [2] and TRAMO model selection [3]. Stock element functions create indi vidual parts of predeﬁned models, they range from the general to the speciﬁc. These functions have names of the form (element) (description) , for example, the function mat arma creates matrices for the ARMA model. Currently there are ﬁ ve types of stock element functions: fun * creates stationary or dynamic matrix update functions, mat * creates matrix structures, ngdist * creates non-Gaussian distributions (or the SSM representation of which), ngfun * creates non-Gaussian distrib ution update functions, and x * creates stock regression v ariables such as trading day v ariables. Note that these functions are created to be class independent whenever possible, in that they do not directly construct or output SSM class objects, although the form of the output may imply the SSM class structures (an exception is that fun * functions return a SSPARAM object). This is primarily due to modular considerations and also that these “utility” functions can then be useful to users who (for whatev er reason) do not wish to use the SSM classes. Miscellaneous helper functions provide additional functionality unsuitable to be implemented in the previous two categories. The most important of which is the function setopt , which speciﬁes the global options or settings structure that is used by almost ev ery data analysis function. One time overrides of individual options are possible for each function call, and each option can be individually speciﬁed, examples of possible options include tolerance and verbosity . Other useful functions include ymarray which generates a year, month array from starting year, starting month and number of months and is useful for functions like x td and x ee which tak e year , month arrays as input, randarma which generates a random AR or MA polynomial, and meancov which calculates the mean and cov ariance of a vector sequence. 7 Data Analysis Examples Many of the data analysis examples are based on the data in the book by Durbin and Koopman [1], to facilitate easy comparison and understanding of the method used here. 7.1 Structural time series models In this example data on road accidents in Great Britain [4] is analyzed using structural time series models following [1]. The purpose the the analysis is to assess the effect of seat belt laws on road accident casualties, with individual monthly ﬁgures for driv ers, front seat passengers and rear seat passengers. The monthly price of petrol and av erage number of kilometers tra veled will be used as re gression v ariables. The data is from January 1969 to December 1984. 16 Peng & Aston 7.1.1 Univariate analysis The dri vers series will be analyzed with univ ariate structural time series model, which consists of local lev el compo- nent µ t , trigonometric seasonal component γ t , regression component (on price of petrol) and intervention component (introduction of seat belt law) β x t . The model equation is y t = µ t + γ t + β x t + ε t , (7) where ε t is the observation noise. The following code e xample constructs this model: lvl = ssmodel(’llm’); seas = ssmodel(’seasonal’, ’trig1’, 12); intv = ssmodel(’intv’, n, ’step’, 170); reg = ssmodel(’reg’, petrol, ’price of petrol’); bstsm = [lvl seas intv reg]; bstsm.name = ’Basic structural time series model’; bstsm.param = [10 0.1 0.001]; The MA TLAB display for object bstsm is bstsm = Basic structural time series model ================================== p = 1 m = 14 r = 12 n = 192 Components (5): ================== [Gaussian noise] p = 1 [local polynomial trend] d = 0 [seasonal] subtype = trig1 s = 12 [intervention] subtype = step tau = 170 [regression] variable = price of petrol H matrix (1, 1) ------------------- SSM T oolbox for Matlab 17 10.0000 Z matrix (1, 14, 192) ---------------------- Stationary part: 1 1 0 1 0 1 0 1 0 1 0 1 DYN DYN Dynamic part: Columns 1 through 9 0 0 0 0 0 0 0 0 0 -2.2733 -2.2792 -2.2822 -2.2939 -2.2924 -2.2968 -2.2655 -2.2626 -2.2655 Columns 10 through 18 0 0 0 0 0 0 0 0 0 -2.2728 -2.2757 -2.2828 -2.2899 -2.2956 -2.3012 -2.3165 -2.3192 -2.3220 ... Columns 181 through 189 1 1 1 1 1 1 1 1 1 -2.1390 -2.1646 -2.1565 -2.1597 -2.1644 -2.1648 -2.1634 -2.1646 -2.1707 Columns 190 through 192 1 1 1 -2.1502 -2.1539 -2.1536 T matrix (14, 14) ------------------- Columns 1 through 9 1 0 0 0 0 0 0 0 0 0 0.8660 0.5000 0 0 0 0 0 0 0 -0.5000 0.8660 0 0 0 0 0 0 0 0 0 0.5000 0.8660 0 0 0 0 0 0 0 -0.8660 0.5000 0 0 0 0 18 Peng & Aston 0 0 0 0 0 0.0000 1 0 0 0 0 0 0 0 -1 0.0000 0 0 0 0 0 0 0 0 0 -0.5000 0.8660 0 0 0 0 0 0 0 -0.8660 -0.5000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Columns 10 through 14 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.8660 0.5000 0 0 0 -0.5000 -0.8660 0 0 0 0 0 -1 0 0 0 0 0 1 0 0 0 0 0 1 R matrix (14, 12) ------------------- 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Q matrix (12, 12) SSM T oolbox for Matlab 19 ------------------- Columns 1 through 9 0.1000 0 0 0 0 0 0 0 0 0 0.0010 0 0 0 0 0 0 0 0 0 0.0010 0 0 0 0 0 0 0 0 0 0.0010 0 0 0 0 0 0 0 0 0 0.0010 0 0 0 0 0 0 0 0 0 0.0010 0 0 0 0 0 0 0 0 0 0.0010 0 0 0 0 0 0 0 0 0 0.0010 0 0 0 0 0 0 0 0 0 0.0010 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Columns 10 through 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0010 0 0 0 0.0010 0 0 0 0.0010 P1 matrix (14, 14) ------------------- Inf 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Inf 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Inf 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Inf 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Inf 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Inf 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Inf 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Inf 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Inf 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Inf 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Inf 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Inf 0 0 20 Peng & Aston 0 0 0 0 0 0 0 0 0 0 0 0 Inf 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Inf Adjacency matrix (3 functions): ----------------------------------------------------------------------------- Function H Z T R Q c a1 P1 Hd Zd Td Rd Qd cd Hng Qng Znl Tnl fun_var/psi2var 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 fun_var/psi2var 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 fun_dupvar/psi2dupvar1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 Parameters (3) ------------------------------------ Name Value epsilon var 10 zeta var 0.1 omega var 0.001 The constituting components of the model and the model matrices are all displayed, note that Z is dynamic with the intervention and regression variables, and matrices c and a1 are omitted if they are zero, also H and Q take on values determined by the parameters set. W ith the model constructed, estimation can proceed with the code: [bstsm logL] = estimate(y, bstsm); [alphahat V] = statesmo(y, bstsm); irr = disturbsmo(y, bstsm); ycom = signal(alphahat, bstsm); ylvl = sum(ycom([1 3 4], :), 1); yseas = ycom(2, :); The estimated model parameters are [0.0037862 0.00026768 1.162e-006] , which can be obtained by dis- playing bstsm.param , the loglikelihood is 175.7790 . State and disturbance smoothing is performed with the estimated model, and the smoothed state is transformed into signal components by signal . Because p = 1 , the out- put ycom is M × n , where M is the number of signal components. The lev el, interv ention and regression are summed as the total data lev el, separating seasonal inﬂuence. Using MA TLAB graphic functions the indi vidual signal components and data are visualized: subplot(3, 1, 1), plot(time, y, ’r:’, ’DisplayName’, ’drivers’), hold all, plot(time, ylvl, ’DisplayName’, ’est. level’), hold off, title(’Level’), ylim([6.875 8]), legend(’show’); subplot(3, 1, 2), plot(time, yseas), title(’Seasonal’), ylim([-0.16 0.28]); subplot(3, 1, 3), plot(time, irr), title(’Irregular’), ylim([-0.15 0.15]); The graphical output is shown in ﬁgure 1. The coefﬁcient for the intervention and regression is deﬁned as part of the state vector in this model, so they can be obtained from the last two elements of the smoothed state vector (due to the order of component concatenation) at the last time point. Coefﬁcient for the intervention is alphahat(13, end) = -0.23773 , and coef ﬁcient for the regression (price of petrol) is alphahat(14, end) = -0.2914 . In this way diagnostics of these coef ﬁcient estimates can also be obtained by the smoothed state variance V . SSM T oolbox for Matlab 21 Figure 1: Driver casualties estimated by basic structural time series model 7.1.2 Bivariate analysis The front seat passenger and rear seat passenger series will be analyzed together using biv ariate structural time series model, with components as before. Speciﬁcally , separate lev el and seasonal components are deﬁned for both series, but the disturbances are assumed to be correlated. T o reduce the number of parameters estimated the seasonal component is assumed to be ﬁxed, so that the total number of parameters is six. W e also include regression on the price of petrol and kilometers traveled, and intervention for only the ﬁrst series, since the seat belt law only effects the front seat passengers. The following is the model construction code: bilvl = ssmodel(’mvllm’, 2); biseas = ssmodel(’mvseasonal’, 2, [], ’trig fixed’, 12); biintv = ssmodel(’mvintv’, 2, n, {’step’ ’null’}, 170); bireg = ssmodel(’mvreg’, 2, [petrol; km]); bistsm = [bilvl biseas biintv bireg]; bistsm.name = ’Bivariate structural time series model’; The model is then estimated with carefully chosen initial values, and state smoothing and signal extraction proceeds as before: [bistsm logL] = estimate(y2, bistsm, [0.0054 0.0086 0.0045 0.00027 0.00024 0.00023]); [alphahat V] = statesmo(y2, bistsm); y2com = signal(alphahat, bistsm); y2lvl = sum(y2com(:,:, [1 3 4]), 3); y2seas = y2com(:,:, 2); When p > 1 the output from signal is of dimension p × n × M , where M is the number of components. The lev el, regression and intervention are treated as one component of data lev el as before, separated from the seasonal component. The components estimated for the two series is plotted: subplot(2, 1, 1), plot(time, y2lvl(1, :), ’DisplayName’, ’est. level’), hold all, 22 Peng & Aston Figure 2: Passenger casualties estimated by biv ariate structural time series model scatter(time, y2(1, :), 8, ’r’, ’s’, ’filled’, ’DisplayName’, ’front seat’), hold off, title(’Front seat passenger level’), xlim([68 85]), ylim([6 7.25]), legend(’show’); subplot(2, 1, 2), plot(time, y2lvl(2, :), ’DisplayName’, ’est. level’), hold all, scatter(time, y2(2, :), 8, ’r’, ’s’, ’filled’, ’DisplayName’, ’rear seat’), hold off, title(’Rear seat passenger level’), xlim([68 85]), ylim([5.375 6.5]), legend(’show’); The graphical output is shown in ﬁgure 2. The intervention coefﬁcient for the front passenger series is obtained by alphahat(25, end) = -0.30025 , the next four elements are the coefﬁcients of the regression of the two series on the price of petrol and kilometers trav eled. 7.2 ARMA models In this example the dif ference of the number of users logged on an internet server [5] is analyzed by ARMA models, and model selection via BIC and missing data analysis are demonstrated. T o select an appropriate ARMA ( p, q ) model for the data various values for p and q are tried, and the BIC of the estimated model for each is recorded, the model with the lowest BIC v alue is chosen. for p = 0 : 5 for q = 0 : 5 arma = ssmodel(’arma’, p, q); [arma logL output] = estimate(y, arma, 0.1); BIC(p+1, q+1) = output.BIC; end end [m i] = min(BIC(:)); arma = estimate(y, ssmodel(’arma’, mod(i-1, 6), floor((i-1)/6)), 0.1); The BIC values obtained for each model is as follo ws: SSM T oolbox for Matlab 23 Figure 3: Forecast using ARMA (1 , 1) model with 50% conﬁdence interval q = 0 q = 1 q = 2 q = 3 q = 4 q = 5 p = 0 6.3999 5.6060 5.3299 5.3601 5.4189 5.3984 p = 1 5.3983 5.2736 5.3195 5.3288 5.3603 5.3985 p = 2 5.3532 5.3199 5.3629 5.3675 5.3970 5.4436 p = 3 5.2765 5.3224 5.3714 5.4166 5.4525 5.4909 p = 4 5.3223 5.3692 5.4142 5.4539 5.4805 5.4915 p = 5 5.3689 5.4124 5.4617 5.5288 5.5364 5.5871 The model with the lo west BIC is ARMA (1 , 1) , second lo west is ARMA (3 , 0) , the former model is chosen for subse- quent analysis. Next missing data is simulated by setting some time points to NaN , model and state estimation can still proceed normally with missing data present. y([6 16 26 36 46 56 66 72 73 74 75 76 86 96]) = NaN; arma = estimate(y, arma, 0.1); yf = signal(kalman(y, arma), arma); Forecasting is equiv alent to treating future data as missing, thus the data set y is appended with as many NaN values as the steps ahead to forecast. Using the previous estimated ARMA (1 , 1) model the Kalman ﬁlter will then effecti vely predict future data points. [a P v F] = kalman([y repmat(NaN, 1, 20)], arma); yf = signal(a(:, 1:end-1), arma); conf50 = 0.675 * realsqrt(F); conf50(1:n) = NaN; plot(yf), hold all, plot([yf+conf50; yf-conf50]’, ’g:’), scatter(1:length(y), y, 10, ’r’, ’s’, ’filled’), hold off, ylim([-15 15]); The plot of the forecast and its conﬁdence interval is sho wn in ﬁgure 3. Note that if the Kalman ﬁlter is replaced with the state smoother , the forecasted v alues will still be the same. 24 Peng & Aston Figure 4: Cubic spline smoothing of motorcycle acceleration data 7.3 Cubic spline smoothing The general state space formulation can also be used to do cubic spline smoothing, by putting the cubic spline into an equiv alent state space form, and accounting for the continuous nature of such smoothing procedures. Here the continuous acceleration data of a simulated motorcycle accident [6] is smoothed by the cubic spline model, which is predeﬁned. spline = estimate(y, ssmodel(’spline’, delta), [1 0.1]); [alphahat V] = statesmo(y, spline); conf95 = squeeze(1.96 * realsqrt(V(1, 1, :)))’; [eps eta epsvar] = disturbsmo(y, spline); The smoothed data and standardized irregular is plotted and sho wn in ﬁgure 4. subplot(2, 1, 1), plot(time, alphahat(1, :), ’b’), hold all, plot(time, [alphahat(1, :) + conf95; alphahat(1, :) - conf95], ’b:’), scatter(time, y, 10, ’r’, ’s’, ’filled’), hold off, title(’Spline and 95% confidence intervals’), ylim([-140 80]), set(gca,’YGrid’,’on’); subplot(2, 1, 2), scatter(time, eps./realsqrt(epsvar), 10, ’r’, ’s’, ’filled’), title(’Standardized irregular’), set(gca,’YGrid’,’on’); It is seen from ﬁgure 4 that the irregular may be heteroscedastic, an easy ad hoc solution is to model the changing variance of the irregular by a continuous version of the local level model. Assume the irregular variance is σ 2 ε h 2 t at time point t and h 1 = 1 , then we model the absolute value of the smoothed irregular abs(eps) with h t as its le vel. The continuous local lev el model with le vel h t needs to be constructed from scratch. contllm = ssmodel(’’, ’continuous local level’, ... 0, 1, 1, 1, ssmat(0, [], true, zeros(1, n), true), ... ’Qd’, {@(X) exp(2 * X) * delta}, {[]}, ssparam({’zeta var’}, ’1/2 log’)); SSM T oolbox for Matlab 25 Figure 5: Cubic spline smoothing of motorcycle acceleration data with heteroscedastic noise contllm = [ssmodel(’Gaussian’) contllm]; alphahat = statesmo(abs(eps), estimate(abs(eps), contllm, [1 0.1])); h2 = (alphahat/alphahat(1)).ˆ2; h2 is then the relati ve magnitude of the noise variance h 2 t at each time point, which can be used to construct a custom dynamic observation noise model as follo ws. hetnoise = ssmodel(’’, ’Heteroscedastic noise’, ... ssmat(0, [], true, zeros(1, n), true), zeros(1, 0), [], [], [], ... ’Hd’, {@(X) exp(2 * X) * h2}, {[]}, ssparam({’epsilon var’}, ’1/2 log’)); hetspline = estimate(y, [hetnoise spline], [1 0.1]); [alphahat V] = statesmo(y, hetspline); conf95 = squeeze(1.96 * realsqrt(V(1, 1, :)))’; [eps eta epsvar] = disturbsmo(y, hetspline); The smoothed data and standardized irregular with heteroscedastic assumption is shown in ﬁgure 5, it is seen that the conﬁdence interval shrank, especially at the start and end of the series, and the irre gular is slightly more uniform. In summary the motorcycle series is ﬁrst modeled by a cubic spline model with Gaussian irregular assumption, then the smoothed irregular magnitude itself is modeled with a local level model. Using the irregular lev el estimated at each time point as the relativ e scale of irregular variance, a ne w h eteroscedastic irregular continuous model is constructed with the estimated level built-in, and plugged into the cubic spline model to obtain new estimates for the motorcycle series. 7.4 Poisson distrib ution error model The road accident casualties and seat belt la w data analyzed in section 7.1 also contains a monthly van dri ver casualties series. Due to the smaller numbers of van driv er casualties the Gaussian assumption is not justiﬁed in this case, 26 Peng & Aston Figure 6: V an driv er casualties and estimated lev el previous methods cannot be applied. Here a poisson distribution is assumed for the data, the mean is exp( θ t ) and the observation distrib ution is p ( y t | θ t ) = exp  θ T t y t − exp( θ t ) − log y t !  The signal θ t in turn is modeled with the structural time series model. The total model is then constructed by con- catenating a poisson distrib ution model with a standard structural time s eries model, the former model will replace the default Gaussian noise model. pbstsm = [ssmodel(’poisson’) ssmodel(’llm’) ... ssmodel(’seasonal’, ’dummy fixed’, 12) ssmodel(’intv’, n, ’step’, 170)]; pbstsm.name = ’Poisson basic STSM’; Model estimation will automatically calculate the Gaussian approximation to the poisson model. Since this is an exponential family distrib ution the data y t also need to be transformed to ˜ y t , which is stored in the output argument output , and used in place of y t for all functions implementing linear Gaussian (Kalman ﬁlter related) algorithms. The following is the code for model and state estimation [pbstsm logL output] = estimate(y, pbstsm, 0.0006, [], ’fmin’, ’bfgs’, ’disp’, ’iter’); [alphahat V] = statesmo(output.ytilde, pbstsm); thetacom = signal(alphahat, pbstsm); Note that the original data y is input to the model estimation routine, which also calculates the transform ytilde . The model estimated then has its Gaussian approximation b uilt-in, and will be treated by the state smoother as a linear Gaussian model, hence the transformed data ytilde needs to be used as input. The signal components thetacom obtained from the smoothed state alphahat is the components of θ t , the mean of y can then be estimated by exp(sum(thetacom, 1)) . The e xponentiated le vel exp( θ t ) with the seasonal component eliminated is compared to the original data in ﬁgure 6. The effect of the seat belt law can be clearly seen. thetalvl = thetacom(1, :) + thetacom(3, :); plot(time, exp(thetalvl), ’DisplayName’, ’est. level’), hold all, plot(time, y, ’r:’, ’DisplayName’, ’data’), hold off, ylim([1 18]), legend(’show’); SSM T oolbox for Matlab 27 Figure 7: t-distribution irregular 7.5 t-distribution models In this example another kind of non-Gaussian models, t-distribution models is used to analyze quarterly demand for gas in UK [7]. A structural time series model with a t-distrib ution hea vy-tailed observ ation noise is constructed similar to the last section, and model estimation and state smoothing performed. tstsm = [ssmodel(’t’) ssmodel(’llt’) ssmodel(’seasonal’, ’dummy’, 4)]; tstsm = estimate(y, tstsm, [0.0018 4 7.7e-10 7.9e-6 0.0033], [], ’fmin’, ’bfgs’, ’disp’, ’iter’); [alpha irr] = fastsmo(y, tstsm); plot(time, irr), title(’Irregular component’), xlim([1959 1988]), ylim([-0.4 0.4]); Since t-distribution is not an exponential family distribution, data transformation is not necessary , and y is used throughout. A plot of the irregular in ﬁgure 7 shows that potential outliers with respect to the Gaussian assumption has been detected by the use of heavy-tailed distrib ution. 7.6 Hillmer -Tiao decomposition In this example seasonal adjustment is performed by the Hillmer-T iao decomposition [2] of airline models. The data (Manufacturing and reproducing magnetic and optical media, U.S. Census Bureau) is ﬁtted with the airline model, then the estimated model is Hillmer-T iao decomposed into an ARIMA components model with trend and seasonal compo- nents. The same is done for the generalized airline model 1 [8], and the seasonal adjustment results are compared. The following are the code to perform seasonal adjustment with the airline model: air = estimate(y, ssmodel(’airline’), 0.1); aircom = ssmhtd(air); ycom = signal(statesmo(y, aircom), aircom); airseas = ycom(2, :); aircom is the decomposed ARIMA components model corresponding to the estimated airline model, and airseas is the seasonal component, which will be subtracted out of the data y to obtain the seasonal adjusted series. ssmhtd automatically decompose ARIMA type models into trend, seasonal and irregular components, plus any extra MA components as permissible. The same seasonal adjustment procedure is then done with the generalized airline model, using parameter estimates from the airline model as initial parameters: 1 Currently renamed Frequency Speciﬁc ARIMA model. 28 Peng & Aston Figure 8: Seasonal adjustment with airline and generalized airline models param0 = air.param([1 2 2 3]); param0(1:3) = -param0(1:3); param0(2:3) = param0(2:3).ˆ(1/12); ga = estimate(y, ssmodel(’genair’, 3, 5, 3), param0); gacom = ssmhtd(ga); ycom = signal(statesmo(y, gacom), gacom); gaseas = ycom(2, :); The code creates a generalized airline 3-5-1 model, Hillmer -T iao decomposition produces the same components as for the airline model since the two models have the same order . From the code it can be seen that the various functions work transparently across different ARIMA type models. Figure 8 shows the comparison between the two seasonal adjustment results. Refer ences [1] J. Durbin and S. J. K oopman. T ime Series Analysis by State Space Methods . Oxford University Press, Ne w Y ork, 2001. [2] S. C. Hillmer and G. C. T iao. An arima-model-based approach to seasonal adjustment. J. Am. Stat. Assoc. , 77:63–70, 1982. [3] V . Gomez and A. Marav all. Automatic modelling methods for univ ariate series. In D. Pena, G. C. T iao, and R. S. Tsay , editors, A Course in T ime Series Analysis . John W iley & Sons, Ne w Y ork, 2001. [4] A. C. Harvey and J. Durbin. The ef fects of seat belt le gislation on british road casualties: A case study in structural time series modelling. J. Roy . Stat. Soc. A , 149:187–227, 1986. SSM T oolbox for Matlab 29 [5] S. Makridakis, S. C. Wheelwright, and R. J. Hyndman. F or ecasting: Methods and Applications, 3r d Ed. John W iley and Sons, Ne w Y ork, 1998. [6] B. W . Silverman. Some aspects of the spline smoothing approach to non-parametric regression curve ﬁtting. J. Roy . Stat. Soc. B , 47:1–52, 1985. [7] S. J. K oopman, A. C. Harve y , J. A. Doornik, and N. Shephard. ST AMP 6.0: Structural T ime Series Analyser , Modeller and Pr edictor . Timberlake Consultants, London, 2000. [8] J.A.D. Aston, D.F . Findley , T .S. McElroy , K.C. Wills, and D.E.K. Martin. Ne w arima models for seasonal time series and their application to seasonal adjustment and forecasting. Resear ch Report Series, US Census Bureau , 2007. A Pr edeﬁned Model Reference The predeﬁned models can be organized into two cate gories, observ ation disturbance models and normal models. The former contains only speciﬁcation of the observation disturbance, and is used primarily for model combination; the latter contains all other models, and can in turn be partitioned into structural time series models, ARIMA type models, and other models. The following is a list of each model, their string code and the arguments required to construct them via the constructor ssmodel . • Observation disturbance models – Gaussian noise: ∼ 2 (’Gaussian’ | ’normal’[, p, cov]) p is the number of variables, def ault 1 . cov is true if they are correlated, default true . – Null noise: ∼ (’null’[, p]) p is the number of variables, def ault 1 . – Poisson err or: ∼ (’poisson’) – Binary error: ∼ (’binary’) – Binomial error: ∼ (’binomial’, k) k is the number of trials, can be a scalar or row v ector . – Negative binomial err or: ∼ (’negbinomial’, k) k is the number of trials, can be a scalar or row v ector . – Exponential error: ∼ (’exp’) – Multinomial error: ∼ (’multinomial’, h, k) h is the number of cells. k is the number of trials, can be a scalar or row v ector . – Exponential family error: ∼ (’expfamily’, b, d2b, id2bdb, c) p ( y | θ ) = exp  y T θ − b ( θ ) + c ( y )  b is the function b ( θ ) . d2b is ¨ b ( θ ) , the second deriv ative of b ( θ ) . 2 Substitute with name of relev ant the class constructor . 30 Peng & Aston id2bdb is ¨ b ( θ ) − 1 ˙ b ( θ ) . c is the function c ( y ) . – Zero-mean stochastic v olatility error: ∼ (’zmsv’) – t-distribution noise: ∼ (’t’[, nu]) nu is the degree of freedom, will be estimated as model parameter if not speciﬁed. – Gaussian mixture noise: ∼ (’mix’) – General error noise: ∼ (’error’) • Structural time series models – Integrated random walk: ∼ (’irw’, d) d is the order of integration. – Local polynomial trend: ∼ (’lpt’, d) d is the order of the polynomial. – Local level model: ∼ (’llm’) – Local level tr end: ∼ (’llt’) – Seasonal components: ∼ (’seasonal’, type, s) type can be ’dummy’ , ’dummy fixed’ , ’h&s’ , ’trig1’ , ’trig2’ or ’trig fixed’ . s is the seasonal period. – Cycle component: ∼ (’cycle’) – Regression components: ∼ (’reg’, x[, varname]) – Dynamic regr ession components: ∼ (’dynreg’, x[, varname]) x is a m × n matrix, m is the number of regression v ariables. varname is the name of the variables. – Intervention components: ∼ (’intv’, n, type, tau) n is the total time duration. type can be ’step’ , ’pulse’ , ’slope’ or ’null’ . tau is the onset time. – Constant components: ∼ (’constant’) – T rading day variables: ∼ (’td6’ | ’td1’, y, m, N) ’td6’ creates a six-variable trading day model. ’td1’ creates a one-variable trading day model. – Length-of-month variables: ∼ (’lom’, y, m, N) – Leap-year variables: ∼ (’ly’, y, m, N) – Easter effect variables: ∼ (’ee’, y, m, N, d) y is the starting year . m is the starting month. N is the total number of months. d is the number of days before Easter . SSM T oolbox for Matlab 31 – Structural time series models: ∼ (’stsm’, lvl, seas, s[, cycle, x]) lvl is ’level’ or ’trend’ . seas is the seasonal type (see seasonal components). s is the seasonal period. cycle is true if there’ s a cycle component in the model, default false . x is explanatory (re gression) variables (see re gression components). – Common levels models: ∼ (’commonlvls’, p, A ast, a ast) y t = " 0 a ∗ # + " I r A ∗ # µ ∗ t + ε t µ ∗ t +1 = µ ∗ t + η ∗ t p is the number of variables (length of y t ). A ast is A ∗ , a ( p − r ) × r matrix. a ast is a ∗ , a ( p − r ) × 1 vector . – Multivariate local le vel models: ∼ (’mvllm’, p[, cov]) – Multivariate local le vel tr end: ∼ (’mvllt’, p[, cov]) – Multivariate seasonal components: ∼ (’mvseasonal’, p, cov, type, s) – Multivariate cycle component: ∼ (’mvcycle’, p[, cov]) p is the number of variables. cov is a logical vector that is true for each correlated disturbances, default all true . Other arguments are the same as the uni variate versions. – Multivariate r egression components: ∼ (’mvreg’, p, x[, dep]) dep is a p × m logical matrix that speciﬁes dependence of data v ariables on regression v ariables. – Multivariate inter vention components: ∼ (’mvintv’, p, n, type, tau) – Multivariate structural time series models: ∼ (’mvstsm’, p, cov, lvl, seas, s[, cycle, x, dep]) cov is a logical vector that speciﬁes the co variance structure of each component in turn. • ARIMA type models – ARMA models: ∼ (’arma’, p, q, mean) – ARIMA models: ∼ (’arima’, p, d, q, mean) – Multiplicative seasonal ARIMA models: ∼ (’sarima’, p, d, q, P, D, Q, s, mean) – Seasonal sum ARMA models: ∼ (’sumarma’, p, q, D, s, mean) – SARIMA with Hillmer -Tiao decomposition: ∼ (’sarimahtd’, p, d, q, P, D, Q, s, gauss) p is the order of AR. d is the order of I. q is the order of MA. P is the order of seasonal AR. D is the order of seasonal I. Q is the order of seasonal MA. s is the seasonal period. mean is true if the model has mean. gauss is true if the irregular component is Gaussian. 32 Peng & Aston – Airline models: ∼ (’airline’[, s]) s is the period, default 12. – Generalized airline models: ∼ (’genair’, nparam, nfreq, freq) nparam is the number of parameters, 3 or 4 . nfreq is the size of the largest subset of frequencies sharing the same parameter . freq is an array containing the members of the smallest subset. – ARIMA component models: ∼ (’arimacom’, d, D, s, phi, theta, ksivar) The arguments match those of the function htd , see its description for details. • Other models – Cubic spline smoothing (continuous time): ∼ (’spline’, delta) delta is the time duration of each data point. – 1/f noise models (approximated by AR): ∼ (’1/f noise’, m) m is the order of the approximating AR process. B Class Refer ence B.1 SSMA T The class SSMAT represents state space matrices, can be stationary or dynamic. SSMAT can be horizontally , v ertically , and block diagonally concatenated with each other, as well as with ordinary matrices. Parenthesis reference can also be used, where two or less indices indicate the stationary part, and three indices indicates a 3-d matrix with time as the third dimension. B.1.1 Subscripted refer ence and assignment .mat RA 3 Stationary part of SSMAT Assignment cannot change size. .mmask R V ariable mask A logical matrix the same size as mat . .n R T ime duration Is equal to size(dvec, 2) . .dmmask R Dynamic mask A logical matrix the same size as mat . .d R Number of dynamic elements Is equal to size(dvec, 1) . .dvec RA Dynamic vector sequence A d × n matrix. Assignment cannot change d . .dvmask R Dynamic vector v ariable mask A d × 1 logical vector . .const R T rue if SSMAT is constant A logical scalar . .sta R T rue if SSMAT is stationary A logical scalar . (i, j) RA Stationary part of SSMAT Assignment cannot change size. (i, j, t) RA SSMAT as a 3-d matrix Assignment cannot change size. B.1.2 Class methods • ssmat SSMAT constructor . SSMAT(m[, mmask]) creates a SSMAT . m can be a 2-d or 3-d matrix, where the third dimension is time. mmask is a 2-d logical matrix masking the v ariable part of matrix, which can be omitted or set to [] for 3 R - reference, A - assignment. SSM T oolbox for Matlab 33 constant SSMAT . SSMAT(m, mmask, dmmask[, dvec, dvmask]) creates a dynamic SSMAT . m is a 2-d matrix forming the stationary part of the SSMAT . mmask is a logical matrix masking the variable part of matrix, which can be set to [] for constant SSMAT . dmmask is a logical matrix masking the dynamic part of SSMAT . dvec is a nnz(dmmask) × n matrix, dvec(:, t) is the values for the dynamic part at time t ( m(dmmask) = dvec(:, t) for time t ), default is a nnz(dmmask) × 1 zero matrix. dvmask is a nnz(dmmask) × 1 logical vector masking the v ariable part of dvec . • isconst ISCONST( · 4 ) returns true if the stationary part of SSMAT is constant. • isdconst ISDCONST( · ) returns true if the dynamic part of SSMAT is constant. • issta ISSTA( · ) returns true if SSMAT is stationary . • size [ ... ] = SIZE( · , ...) is equiv alent to [ ... ] = SIZE( · .mat, ...) . • getmat GETMAT( · ) returns SSMAT as a matrix if stationary , or as a cell array of matrices if dynamic. • getdvec GETDVEC( · , mmask) returns the dynamic elements speciﬁed in matrix mask mmask as a vector sequence. • getn GETN( · ) returns the time duration of SSMAT . • setmat SETMAT( · , vec) updates the variable stationary matrix, vec must ha ve size nnz(mmask) × 1 . • setdvec SETDVEC( · , dsubvec) updates the variable dynamic vector sequence, dsubvec must hav e nnz(dvmask) rows. • plus SSMAT addition. B.2 SSDIST The class SSDIST represents state space distributions, a set of non-Gaussian distributions that governs a disturbance vector (which can also ha ve Gaussian elements), and is a child class of SSMAT . SSDIST are constrained to be square matrices, hence can only be block diagonally concatenated, any other types of concatenation will ignore the non- Gaussian distribution part. Predeﬁned non-Gaussian distributions can be constructed via the SSDIST constructor . 4 An instance of the class in question. 34 Peng & Aston B.2.1 Subscripted refer ence and assignment .nd R Number of distributions A scalar . .type R T ype of each distribution A 0-1 row v ector . .matf R Approximation functions A cell array of functions. .logpf R Log probability functions A cell array of functions. .diagmask R Non-Gaussian masks A size(mat, 1) × nd logical matrix. .dmask R V ariable mask for the distrib utions A 1 × nd logical vector . .const R T rue if SSDIST is constant A logical scalar . B.2.2 Class methods • ssdist SSDIST constructor . SSDIST(’’, type[, matf, logpf]) creates a single univ ariate non-Gaussian distribution. type is 0 for exponential family distributions and 1 for additi ve noise distributions. matf is the function that generates the approximated Gaussian variance matrix given observations and signal or disturbances. logpf is the function that calculates the log probability of observations given observation and signal or disturbances. If the last two arguments are omitted, then the SSDIST is assumed to be variable. SSDIST with multiple non-Gaussian distributions can be formed by constructing a SSDIST for each and then block diagonally concatenating them. • isdistconst ISDISTCONST( · ) returns true if the non-Gaussian part of SSDIST is constant. • setdist SETDIST( · , distf) updates the non-Gaussian distributions. distf is a cell matrix of functions. • setgauss SETGAUSS( · , eps) or [ · ytilde] = SETGAUSS( · , y, theta) calculates the approximating Gaussian cov ariance, the ﬁrst form is used if there are no exponential f amily distributions, and the data y need to be transformed into ytilde in the second form if there are. • logprobrat LOGPROBRAT( · , N, eps) or LOGPROBRAT( · , N, y, theta, eps) calculates the log probability ratio between the original non-Gaussian distribution and the approximating Gaussian distribution of the observations, it is used when calculating the non-Gaussian loglikelihood by importance sampling. The ﬁrst form is used if there are no exponential family distributions. N is the number of importance samples. B.2.3 Predeﬁned non-Gaussian distrib utions • Poisson distrib ution: ∼ (’poisson’) • Binary distribution: ∼ (’binary’) • Binomial distribution: ∼ (’binomial’, k) • Negative binomial distrib ution: ∼ (’negbinomial’, k) k is the number of trials, can be a scalar or row v ector . SSM T oolbox for Matlab 35 • Exponential distribution: ∼ (’exp’) • Multinomial distribution: ∼ (’multinomial’, h, k) h is the number of cells. k is the number of trials, can be a scalar or row v ector . • General exponential family distribution: ∼ (’expfamily’, b, d2b, id2bdb, c) p ( y | θ ) = exp  y T θ − b ( θ ) + c ( y )  b is the function b ( θ ) . d2b is ¨ b ( θ ) , the second deriv ative of b ( θ ) . id2bdb is ¨ b ( θ ) − 1 ˙ b ( θ ) . c is the function c ( y ) . B.3 SSFUNC The class SSFUNC represents state space functions, a set of nonlinear functions that describe state vector transfor- mations, and is a child class of SSMAT . Since linear transformations in SSM are represented by left multiplication of matrices, it is only possible for nonlinear functions (and hence it’ s matrix approximation) to be either block diagonally or horizontally concatenated. The former concatenation method produces multiple resulting output vectors, whereas the latter produces the sum of the output vectors of all functions as output. V ertical concatenation will ignore the nonlinear part of SSFUNC . B.3.1 Subscripted refer ence and assignment .nf R Number of functions A scalar . .f R The nonlinear functions A cell array of functions. .df R Deriv ative of each function A cell array of functions. .horzmask R Nonlinear function output masks A size(mat, 1) × nf logical matrix. .vertmask R Nonlinear function input masks A size(mat, 2) × nf logical matrix. .fmask R V ariable mask for the functions A 1 × nf logical vector . .const R T rue if SSFUNC is constant A logical scalar . B.3.2 Class methods • ssfunc SSFUNC constructor . SSFUNC(’’, p, m[, f, df]) creates a single nonlinear function. p is the input vector dimension and m is the output vector dimension. Hence the approximating linear matrix will be p × m . f and df is the nonlinear function and its deriv ativ e. f maps m × 1 vectors and time t to p × 1 vectors, and df maps m × 1 vectors and time t to p × m matrices. If f and df are not both provided, the SSFUNC is assumed to be variable. SSFUNC with multiple functions can be constructed by deﬁning indi vidual SSFUNC for each function, and concatenating them. • isfconst ISFCONST( · ) returns true if the nonlinear part of SSFUNC is constant. 36 Peng & Aston • getfunc GETFUNC( · , x, t) returns f(x, t) , the transform of x at time t . • getvertmask GETVERTMASK( · ) returns vertmask . • setfunc SETFUNC( · , funcf) updates the nonlinear functions. funcf is a cell matrix of functions. • setlinear [ · c] = SETLINEAR( · , alpha) calculates the approximating linear matrix, c is an additive constant in the approximation. B.4 SSP ARAM The class SSPARAM represents state space model parameters. SSPARAM can be seen as a row vector of parameter values, hence the y are combined by horizontal concatenation with a slight dif ference, as described below . B.4.1 Subscripted refer ence and assignment .w R Number of parameters A scalar . .name R Parameter names A cell array of strings. .value RA T ransformed parameter values A 1 × w vector . Assignment cannot change size. B.4.2 Class methods • ssparam SSPARAM constructor . SSPARAM(w[, transform]) creates a SSPARAM with w parameters. SSPARAM(name[, transform]) creates a SSPARAM with parameter names name , which is a cell array of strings. transform is a cell array of strings describing transforms used for each parameter . SSPARAM(name, group, transform) creates a SSPARAM with parameter names speciﬁed by name . group is a row vector that speciﬁes the number of parameters included in each parameter subgroup, and transform is of the same length which speciﬁes transforms for each group. • get GET( · ) returns the untransformed parameter values. • set SET( · , value) sets the untransformed parameter values. • remove REMOVE( · , mask) removes the parameters speciﬁed by mask . • horzcat [ · pmask] = HORZCAT( · 1 , · 2 , ...) combines multiple SSPARAM into one, and optionally generates a cell array of parameter masks pmask , masking each original parameter sets with respect to the ne w combined parameter set. SSM T oolbox for Matlab 37 B.5 SSMODEL The class SSMODEL represents state space models and embeds all the previous classes. SSMODEL can be combined by horizontal or block diagonal concatenation. The former is the basis for additiv e model combination, where the observation is the sum of individual signals. The latter is an independent combination in which the observ ation is the vertical concatenation of the indi vidual signals, which obviously needs to be further combined or changed to be useful. Only SSMODEL objects are used directly in data analysis functions, knowledge of embedded objects are needed only when deﬁning custom models. For a list of predeﬁned models see appendix A. B.5.1 Subscripted refer ence and assignment .name RA Model name A string. .Hinfo R Observation disturbance info A structure. .info R Model component info A cell array of structures. .p R Model observ ation dimension A scalar . .m R Model state dimension A scalar . .r R Model state disturbance dimension A scalar . .n R Model time duration A scalar . V alue is 1 for stationary models. .H R Observ ation disturbance A SSMAT or SSDIST . .Z R State to observ ation transform A SSMAT or SSFUNC . .T R State update transform A SSMAT or SSFUNC . .R R Disturbance to state transform A SSMAT . .Q R State disturbance A SSMAT or SSDIST . .c R State update constant A SSMAT . .a1 RA Initial state mean A stationary SSMAT , can assign matrices. .P1 RA Initial state variance A stationary SSMAT , can assign matrices. .sta R T rue for stationary SSMODEL A logical scalar . .linear R T rue for linear SSMODEL A logical scalar . .gauss R T rue for Gaussian SSMODEL A logical scalar . .w R Number of parameters A scalar . .paramname R Parameter names A cell array of strings. .param RA Parameter v alues A 1 × w v ector . Assignment cannot change w . .psi RA T ransformed parameter values A 1 × w vector . Assignment cannot change w . .q R Number of initial dif fuse elements A scalar . B.5.2 Class Methods • ssmodel SSMODEL constructor . SSMODEL(’’, info, H, Z, T, R, Q) creates a constant state space model with complete diffuse ini- tialization and c = 0 . SSMODEL(’’, info, H, Z, T, R, Q, A, func, grad, psi[, pmask, P1, a1, c]) create a state space model. info is a string or a free form structure with ﬁeld ’type’ describing the model component. H is a SSMAT or SSDIST . Z is a SSMAT or SSFUNC . T is a SSMAT or SSFUNC . R is a SSMAT . Q is a SSMAT or SSDIST . A is a cell array of strings or a single string, each corresponding to func , a cell array of functions or a single function. Each string of A is a concatenation of some of ’H’ , ’Z’ , ’T’ , ’R’ , ’Q’ , ’c’ , ’a1’ , ’P1’ , ’Hd’ , ’Zd’ , ’Td’ , ’Rd’ , ’Qd’ , ’cd’ , ’Hng’ , ’Qng’ , ’Znl’ , ’Tnl’ , 38 Peng & Aston representing each part of each element of SSMODEL . For example, ’Qng’ speciﬁes the non-Gaussian part of Q , and ’Zd’ speciﬁes the dynamic part of Z . A complicated example of a string in A is [’H’ repmat(’ng’, 1, gauss) repmat(’T’, 1, p+P>0) ’RQP1’] , which is the adjacency string for a possibly non-Gaussian ARIMA components model. func is a cell array of functions that updates model matrices. If multiple parts of the model are updated the output for each part must be ordered as for the adjacenc y matrix, but parts not updated can be skipped, for example, a function that updates H , T and Q must output [Hvec Tvec Qvec] in this order . grad is the deri vati ve of func , each of which can be [] if not differentiable. For the example function abov e the gradient function would output [Hvec Tvec Qvec Hgrad Tgrad Qgrad] in order . psi is a SSPARAM . pmask is a cell array of parameter masks for the corresponding functions, can be omitted or set to [] if there’ s only one update function (which presumably uses all parameters). P1 and a1 are both stationary SSMAT . c is a SSMAT . All arguments expecting SSMAT can also take numeric matrices. • isgauss ISGAUSS( · ) is true for Gaussian SSMODEL . • islinear ISLINEAR( · ) is true for linear SSMODEL . • issta ISSTA( · ) is true for stationary SSMODEL . Note that all SSFUNC are implicitly assumed to be dynamic, but are treated as stationary by this and similar functions. • set SET( · , A, M, func, grad, psi, pmask) changes one speciﬁc element of SSMODEL . A is an adja- cency string, specifying which elements the function or functions func updates. A is a adjacency string as described earlier , corresponding to func . M should be the ne w v alue for one of H , Z , T , R , Q , c , a1 or P1 , and A is constrained to contain only references to that element. grad is the same type as func and is its gradient, set individual functions to 0 if gradient does not exist or needed. psi is a SSPARAM that stores all the parameters used by functions in func , and pmask is the parameter mask for each function. • setparam SETPARAM( · , psi, transformed) sets the parameters for SSMODEL . If transformed is true, psi is the new v alue for the transformed parameters, else psi is the new v alue for the untransformed parameters. C Function Refer ence C.1 User -deﬁned update functions This section details the format of update functions for various parts of classes SSMAT , SSDIST and SSFUNC . • Stationary part vec = func(psi) updates the stationary part of SSMAT . vec must be a column vector such that mat(mmask) = vec would correctly update the stationary matrix mat giv en parameters psi . SSM T oolbox for Matlab 39 • Dynamic part dsubvec = func(psi) updates the dynamic part of SSMAT . dsubvec must be a matrix such that dvec(dvmask, 1:size(dsubvec, 2)) = dsubvec would correctly update the dynamic vector sequence dvec given parameters psi . • Non-Gaussian part distf = func(psi) updates the non-Gaussian part of SSDIST . distf must be a cell matrix of function handles such that [matfdmask] = distf:, 1 and [logpfdmask] = distf:, 2 would correctly update the distributions gi ven parameters psi . • Nonlinear part funcf = func(psi) updates the nonlinear part of SSFUNC . funcf must be a cell matrix of function handles such that [ffmask] = funcf:, 1 and [dffmask] = funcf:, 2 would correctly update the nonlinear functions giv en parameters psi . Note that an update function can return more than one of the four kind of output, thus updating multiple part of a SSMODEL , b ut the order of the output arguments must follow the con vention: stationary part of H , Z , T , R , Q , c , a1 , P1 , dynamic part of H , Z , T , R , Q , c , non-Gaussian part of H , Q , nonlinear part of Z , T , with possible omissions. C.2 Data analysis functions Most functions in this section accepts analysis settings options, speciﬁed as option name and option value pairs (e.g. (’disp’, ’off’) ). These groups of arguments are speciﬁed at the end of each function that accepts them, and are represented by opt in this section. • batchsmo [alphahat epshat etahat] = BATCHSMO(y, model[, opt]) performs batch smoothing of mul- tiple data sets. y is the data of dimension p × n × N , where n is the data length and N is the number of data sets, there must be no missing values. model is a SSMODEL . The output is respectively the smoothed state, smoothed observation disturbance and smoothed state disturbance, each of dimensions m × n × N , p × n × N and r × n × N . This is equiv alent to doing fastsmo on each data set. • disturbsmo [epshat etahat epsvarhat etavarhat] = DISTURBSMO(y, model[, opt]) performs dis- turbance smoothing. y is the data of dimension p × n , and model is a SSMODEL . The output is respectiv ely the smoothed observation disturbance ( p × n ), smoothed state disturbance ( r × n ), smoothed observation dis- turbance variance ( p × p × n or 1 × n if p = 1 ) and smoothed state disturbance variance ( r × r × n or 1 × n if r = 1 ). • estimate [model logL output] = ESTIMATE(y, model[, param0, alpha0, opt]) estimates the pa- rameters of model starting from the initial parameter value param0 . y is the data of dimension p × n , and model is a SSMODEL . param0 can be empty if the current parameter v alues of model is used as initial v alue, and a scalar param0 sets all parameters to the same value. Alternati vely param0 can be a logical row vec- tor specifying which parameters to estimate, or a 2 × w matrix with the ﬁrst row as initial value and second row as estimated parameter mask. The initial state sequence estimate alpha0 is needed only when model is non-Gaussian or nonlinear . output is a structure that contains optimization routine information, approximated observation sequence ˜ y if non-Gaussian or nonlinear , and the AIC and BIC of the output model . 40 Peng & Aston • fastsmo [alphahat epshat etahat] = fastsmo(y, model[, opt]) performs fast smoothing. y is the data of dimension p × n , and model is a SSMODEL . The output is respectively the smoothed state ( m × n ), smoothed observation disturbance ( p × n ), and smoothed state disturbance ( r × n ). • gauss [model ytilde] = GAUSS(y, model[, alpha0, opt]) calculates the Gaussian approximation. y is the data of dimension p × n , and model is a SSMODEL . alpha0 is the initial state sequence estimate and can be empty or omitted. • kalman [a P v F] = KALMAN(y, model[, opt]) performs Kalman ﬁltering. y is the data of dimension p × n , and model is a SSMODEL . The output is respectiv ely the ﬁltered state ( m × n+1 ), ﬁltered state variance ( m × m × n+1 , or 1 × n+1 if m = 1 ), one-step prediction error (innovation) ( p × n ), one-step prediction v ariance ( p × p × n , or 1 × n if p = 1 ). • linear [model ytilde] = LINEAR(y, model[, alpha0, opt]) calculates the linear approximation. y is the data of dimension p × n , and model is a SSMODEL . alpha0 is the initial state sequence estimate and can be empty or omitted. • loglik LOGLIK(y, model[, ytilde, opt]) returns the log likelihood of model giv en y . y is the data of dimension p × n , and model is a SSMODEL . ytilde is the approximating observation ˜ y and is needed for some non-Gaussian or nonlinear models. • sample [y alpha eps eta] = SAMPLE(model, n[, N]) generates observation samples from model . model is a SSMODEL , n speciﬁes the sampling data length, and N speciﬁes how many sets of data to generate. y is the sampled data of dimension p × n × N , alpha , eps , eta are respectiv ely the corresponding sampled state ( m × n × N ), observation disturbance ( p × n × N ), and state disturbance ( r × n × N ). • signal SIGNAL(alpha, model) generates the signal for each component according to the state sequence and model speciﬁcation. alpha is the state of dimension m × n , and model is a SSMODEL . The output is a cell array of data each with dimension p × n , or a M × n matrix where M is the number of components if p = 1 . • simsmo [alphatilde epstilde etatilde] = SIMSMO(y, model, N[, antithetic, opt]) gen- erates observation samples from model conditional on data y . y is the data of dimension p × n , and model is a SSMODEL . antithetic should be set to 1 if antithetic v ariables are used. The output is respectiv ely the sampled state sequence ( m × n × N ), sampled observ ation disturbance ( p × n × N ), and sampled state disturbance ( r × n × N ). • statesmo [alphahat V] = STATESMO(y, model[, opt]) performs state smoothing. y is the data of dimen- sion p × n , and model is a SSMODEL . The output is respectiv ely the smoothed state ( m × n ), smoothed state variance ( m × m × n , or 1 × n if p = 1 ). If only the ﬁrst output argument is speciﬁed, fast state smoothing is automatically performed instead. SSM T oolbox for Matlab 41 • arimaselect [p d q mean] = ARIMASELECT(y) or [p d q P D Q mean] = ARIMASELECT(y, s) performs TRAMO model selection. y is the data of dimension p × n , and s is the seasonal period. The ﬁrst form selects among ARIMA models with or without a mean. The second form selects among SARIMA models with or without a mean. • armadegree [p q] = ARMADEGREE(y[, mean, mr]) or [p q P Q] = ARMADEGREE(y, s[, mean, mr, ms]) determines the degrees of ARMA or SARMA from observ ation y . y is the data of dimension p × n , and s is the seasonal period. mean is true if models with mean are used for the selection. mr and ms is the largest degree to consider for regular ARMA and seasonal ARMA respectiv ely . • diffdegree [d mean] = DIFFDEGREE(y[, ub, tsig]) or [d D mean] = DIFFDEGREE(y, s[, ub, tsig]) performs unit root detection. The ﬁrst form only considers regular dif ferencing. y is the data of dimension p × n , s is the seasonal period, ub is the threshold for unit roots, default is [0.97 0.88] , and tsig is the t-value threshold for mean detection, def ault is 1.5 . • htd [theta ksivar] = HTD(d, D, s, Phi, Theta[, etavar]) performs Hillmer-T iao decompo- sition. d is the order of regular differencing, D is the order of seasonal differencing, s is the seasonal period, Phi is a cell array of increasing order polynomials representing autoregressiv e factors for each component, Theta is an increasing order polynomial representing the moving average factor , etavar is the disturbance variance. The output theta is a cell array of decomposed moving av erage factors for each component, and ksivar is the corresponding disturbance v ariance for each component (including the irregular v ariance at the end). • ssmhtd SSMHTD(model) performs Hillmer-T iao decomposition on SSMODEL model and returns an ARIMA com- ponent model. This is a wrapper for the function htd . • loglevel LOGLEVEL(y, s) or LOGLEVEL(y, p, q) or LOGLEVEL(y, p, d, q[, P, D, Q, s]) determines the log le vel speciﬁcation with respect to various models. y is the observation data, s is the seasonal period. The output is 1 if no logs are needed, else 0 . The ﬁrst form uses the airline model with speciﬁed seasonal period, second form uses ARMA models, and the third form uses SARIMA models. • oosforecast [yf err SS] = OOSFORECAST(y, model, n1, h) performs out-of-sample forecast. y is the data of dimension p × n , model is a SSMODEL , n1 is the number of time points to exclude at the end, and h is the number of steps ahead to forecast, which can be an array . The output yf is the forecast obtained, err is the forecast error , and SS is the forecast error cumulati ve sum of squares. 42 Peng & Aston C.3 Stock element functions • fun arma [fun grad psi] = FUN ARMA(p, q) creates matrix update functions for ARMA models. • fun cycle [fun grad psi] = FUN CYCLE() creates matrix update functions for cycle components. • fun dupvar [fun grad psi] = FUN DUPVAR(p, cov, d[, name]) creates matrix update functions for dupli- cated variances. • fun genair [fun grad psi] = FUN GENAIR(nparam, nfreq, freq) creates matrix update functions for gen- eralized airline models. • fun homovar [fun grad psi] = FUN HOMOVAR(p, cov, q[, name]) creates matrix update functions for ho- mogeneous variances. • fun interlvar [fun grad psi] = FUN INTERLVAR(p, q, cov[, name]) creates matrix update functions for q- interleav ed variances. • fun mvcycle [fun grad psi] = FUN MVCYCLE(p) creates matrix update functions for multiv ariate cycle component. • fun oneoverf [fun grad psi] = FUN ONEOVERF(m) creates matrix update functions for 1/f noise. • fun sarimahtd [fun grad psi] = FUN SARIMAHTD(p, d, q, P, D, Q, s, gauss) creates matrix update func- tions for SARIMA with Hillmer-T iao decomposition and non-Gaussian irregular . • fun sarma [fun grad psi] = FUN SARMA(p, q, P, Q, s) creates matrix update functions for seasonal ARMA models. • fun spline [fun grad psi] = FUN SPLINE(delta) creates matrix update functions for the cubic spline models. • fun var [fun grad psi] = FUN VAR(p, cov[, name]) creates matrix update functions for variances. • fun wvar [fun grad psi] = FUN WVAR(p, s[, name]) creates matrix update functions for W structure vari- ances. • mat arima [Z T R P1 Tmmask Rmmask P1mmask] = MAT ARIMA(p, d, q, mean) creates base matrices for ARIMA models. SSM T oolbox for Matlab 43 • mat arma [Z T R P1 Tmmask Rmmask P1mmask] = MAT ARMA(p, q, mean) or [T R P1 Tmmask Rmmask P1mmask] = MAT ARMA(p, q, mean) creates base matrices for ARMA models. • mat commonlvls [Z T R a1 P1] = MAT COMMONLVLS(p, A, a) creates base matrices for common levels models. • mat dummy [Z T R] = MAT DUMMY(s[, fixed]) creates base matrices for dummy seasonal component. • mat dupvar [m mmask] = MAT DUPVAR(p, cov, d) creates base matrices for duplicated variances. • mat homovar [H Q Hmmask Qmmask] = MAT HOMOVAR(p, cov, q) creates base matrices for homogeneous vari- ances. • mat hs [Z T R] = MAT HS(s) creates base matrices for Harrison and Stev ens seasonal component. • mat interlvar [m mmask] = MAT INTERLVAR(p, q, cov) creates base matrices for q-interleaved v ariances. • mat lpt [Z T R] = MAT LPT(d, stochastic) creates base matrices for local polynomial trend or integrated random walk models. • mat mvcycle [Z T Tmmask R] = MAT MVCYCLE(p) creates base matrices for multiv ariate cycle components. • mat mvdummy [Z T R] = MAT MVDUMMY(p, s[, fixed]) creates base matrices for multiv ariate dummy seasonal components. • mat mvhs [Z T R] = MAT MVHS(p, s) creates base matrices for multiv ariate Harrison and Stevens seasonal com- ponents. • mat mvintv [Z Zdmmask Zdvec T R] = MAT MVINTV(p, n, type, tau) creates base matrices for multiv ari- ate intervention components. • mat mvllt [Z T R] = MAT MVLLT(p) creates base matrices for multiv ariate local le vel trend models. • mat mvreg [Z Zdmmask Zdvec T R] = MAT MVREG(p, x, dep) creates base matrices for multiv ariate regres- sion components. 44 Peng & Aston • mat mvtrig [Z T R] = MAT MVTRIG(p, s[, fixed]) creates base matrices for multiv ariate trigonometric sea- sonal components. • mat reg [Z Zdmmask Zdvec T R] = MAT REG(x[, dyn]) creates base matrices for regression components. • mat sarima [Z T R P1 Tmmask Rmmask P1mmask] = MAT SARIMA(p, d, q, P, D, Q, s, mean) or [Z T R P1 Rmmask P1mmask] = MAT SARIMA(p, d, q, P, D, Q, s, mean) creates base ma- trices for SARIMA models. • mat sarimahtd [Z T R Qmat P1 Tmmask Rmmask Qmmask P1mmask] = MAT SARIMAHTD(p, d, q, P, D, Q, s) creates base matrices for SARIMA with Hillmer-T iao decomposition. • mat spline [Z T Tdmmask Tdvec R] = MAT (delta) creates base matrices for the cubic spline model. • mat sumarma [Z T R P1 Tmmask Rmmask P1mmask] = MAT SUMARMA(p, q, D, s, mean) creates base ma- trices for sum integrated ARMA models. • mat trig [Z T R] = MAT TRIG(s[, fixed]) creates base matrices for trigonometric seasonal components. • mat var [m mmask] = MAT VAR(p, cov) creates base matrices for variances. • mat wvar [m mmask] = MAT WVAR(p, s) creates base matrices for W structure variances. • ngdist binomial [matf logpf] = NGDIST BINOMIAL(k) creates base functions for binomial distrib utions. • ngdist expfamily [matf logpf] = NGDIST EXPFAMILY(b, d2b, id2bdb, c) creates base functions for exponen- tial family distrib utions. • ngdist multinomial [matf logpf] = NGDIST MULTINOMIAL(h, k) creates base functions for multinomial distributions. • ngdist negbinomial [matf logpf] = NGDIST NEGBINOMIAL(k) creates base functions for negati ve binomial distrib utions. • ngfun t [fun psi] = NGFUN T([nu]) creates update functions for the t-distrib ution. • psi2err distf = PSI2ERR(X) is the update function for the general error distribution. SSM T oolbox for Matlab 45 • psi2mix distf = PSI2MIX(X) is the update function for the Gaussian mixture distribution. • psi2zmsv distf = PSI2ZMSV(X) is the update function for the zero-mean stochastic volatility model disturbance distribution. • x ee X EE(Y, M, d) generates Easter effect v ariable. • x intv X INTV(n, type, tau) generates intervention variables. • x ly X LY(Y, M) generates leap-year variable. • x td X TD(Y, M, td6) generates trading day variables. C.4 Miscellaneous helper functions • diaginprod y = DIAGINPROD(A, x1[, x2]) calculates the inner product of a vector sequence. A is a m × m square matrix. x1 and x2 are m × n matrices representing a set of n column vectors each of size m × 1 . The ( i, j ) th element of y is equal to (x1(:, j) − x2(:, i)) T A(x1(:, j) − x2(:, i)) . • f alpha2arma [xAR xMA] = F ALPHA2ARMA(n, m, alpha, sigma2) generates 1/f noise by AR and MA approx- imations. n is length of series to generate, m is number of terms to use in the approximation, alpha is the exponent of in verse frequency of the 1/f noise, and sigma2 the variance of the 1/f noise. The output is the 1/f noise generated by AR and MA approximation respectiv ely . • meancov [m C] = MEANCOV(x) calculates the mean and cov ariance of each vector set in a sequence. x is a m × n × N matrix representing n vector sets each of size N containing m × 1 vectors. m is a m × n matrix, the mean of each of the n vector sets. C is a m × m × n matrix, the cov ariance of each of the n vector sets. • randarma RANDARMA(n, r) generates random ARMA polynomials. n is the polynomial degree and r is root magnitude range. • setopt SETOPT(optname1, optvalue1, optname2, optvalue2, ...) sets the global analysis settings and returns a structure containing the ne w settings, thus a call with no arguments just returns the current settings. The av ailable options are: ’disp’ can be set to ’off’ , ’notify’ , ’final’ or ’iter’ . tol is the tolerance. fmin is the minimization algorithm to use ( ’bfgs’ or ’simplex’ ). maxiter is the maximum number of iterations. nsamp is the number of samples to use in simulation. Note that the same form of argument sequence can be used at the end of most data analysis functions to specify one time ov errides. 46 Peng & Aston • ymarray [Y M] = YMARRAY(y, m, N) generates an array of years and months from speciﬁed starting year y , start- ing month m , and number of months N .

The SSM Toolbox for Matlab

Original Paper

Comments & Academic Discussion

Leave a Comment

Original Paper

Related Papers

Comments & Academic Discussion

Leave a Comment