Science, Faculty of
Statistics, Department of
DSpace
UBCV
Mazzi, Sonia V.T.
Doctor of Philosophy - PhD
University of British Columbia
Nonparametric regression methods are devised in order to obtain a smooth fit of a regression curve without having to specify a parametric model. This, in principle, provides a flexible approach to curve estimation. It is shown how splines are not such a flexible tool for smoothing. It is proposed to use the modern approach of state space modeling instead, to solve a variety of smoothing problems. This model based approach alleviates estimation and computational problems as well as provides insightful solutions to smoothing problems.
5017486 bytes
application/pdf
T H E M O D E L B A S E D A P P R O A C H T O S M O O T H I N G by SONIA V . T . M A Z Z I L i e , Universidad Nacional de Cordoba, Argentina, 1989 M . S c , The University of British Columbia, 1992 A T H E S I S S U B M I T T E D IN P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F D O C T O R O F P H I L O S O P H Y in T H E F A C U L T Y O F G R A D U A T E S T U D I E S D E P A R T M E N T O F S T A T I S T I C S We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A January 1997 © S o n i a V . T . Mazzi , 1997 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of €Tr\ Tl'STt'CS The University of British Columbia Vancouver, Canada Date Sartygr^j 39 -No, DE-6 (2/88) The model based approach to smoothing Abstract Nonparametric regression methods are devised in order to obtain a smooth fit of a regression curve without having to specify a parametric model. This, in principle, provides a flexible approach to curve estimation. It is shown how splines are not such a flexible tool for smoothing. It is proposed to use the modern approach of state space modeling instead, to solve a variety of smoothing problems. This model based approach alleviates estimation and computational problems as well as provides insightful solutions to smoothing problems. n Contents Abstract i i Table of Contents i i i List of Tables v i i i List of Figures ix Acknowledgements x 1 Introduction 1 1.1 Motivation 1 1.2 Contributions 6 1.3 Outline 8 2 The State space model 10 2.1 Definition of the state space model 14 iii 2.2 Some examples of state space models 18 2.2.1 The parametric regression model as a S S M 18 2.2.2 Structural models 19 2.2.3 A R I M A models as state space models 24 2.3 Principles and tools of prediction and estimation for the S S M 28 2.3.1 M i n i m u m mean square error prediction 30 2.3.2 Combining predictors and iterated prediction 32 2.3.3 The prediction and Kalman update formulas 32 2.3.4 The Kalman filter 34 2.3.5 The smoothing filter 35 2.4 The diffuse Kalman filter 39 2.4.1 The diffuse likelihood 40 2.4.2 Diffuse prediction 41 2.4.3 Diffuse smoothing 42 2.5 State space models which have applied since time immemorial 43 2.5.1 Prediction and likelihood evaluation for non-stationary models . . . 44 2.6 Cross-validation and related methods of estimation 46 2.6.1 Cross-validation for state space models 47 iv 2.6.2 Generalised cross-validation and state space models 48 2.6.3 Generalised cross-validation and diffuse state space models 49 2.7 Summary 50 3 Continuous state space models 51 3.1 Definition of continuous state space models 53 3.2 Discrete version of a continuous state space model 53 3.3 Estimation using the continuous state space model 54 3.3.1 The state transition matrix 54 3.3.2 Variance matrix of the state error 55 3.3.3 Initialization of the Diffuse Kalman filter when using a continuous state space model 63 3.4 Smoothing L-Splines 67 3.4.1 Definition of L-splines 68 3.4.2 A stochastic model for L-spline smoothing 69 3.5 Some continuous structural models 76 3.6 Summary 82 4 Smoothing polynomial splines 84 v 4.1 Definition of polynomial smoothing splines 86 4.2 A basis for the space of natural splines 89 4.2.1 Divided differences 89 4.2.2 The Demmler-Reinsch basis 91 4.3 Estimation of the smoothing parameter 95 4.4 A Bayesian approach to spline smoothing 96 4.4.1 Some Gaussian processes 96 4.4.2 Wahba's approach 99 4.4.3 A state space form of Wahba's model 101 4.5 Summary 104 A closer look at polynomial splines 105 5.1 The signal plus noise model 106 5.2 Penalised least squares 107 5.3 Smoothing splines, penalized least squares, signal extraction and state space models 108 5.4 Smoothing splines and A R I M A models I l l 5.5 Related topics 114 5.5.1 Smoothing in the presence of heteroskedasticity 114 vi 5.5.2 Repeated measurements 115 5.5.3 Semiparametric regression 115 5.5.4 Change point problems or intervention analysis 116 5.6 Summary 116 Epilogue 117 References 118 A Stochastic differential equations 124 A . l Definition of the stochastic integral 124 A.2 Stochastic differentials 125 A.3 Stochastic differential equations, existence and uniqueness of solutions . . . 126 A.4 Linear stochastic differential equations 127 B Kronecker products 129 C Jordan canonical form of a matrix 131 D The matrix exponential 132 E Data sets 134 vii List of Tables 2.1 Different trend models 23 3.1 Estimated parameters of chlorine data model for two different modeling situations 73 3.2 Estimated parameters for the L-spline and polynomial spline of order 7 models for melanoma data 79 vm List of Figures 1.1 A good data analysis paradigm 4 3.1 The time immemorial assumption 63 3.2 Predicted signals with 95% confidence bands of chlorine data using two different models 74 3.3 Predicted innovation residuals with 95% confidence bands for chlorine data using two different models 75 3.4 Predicted signals for Case I and II, superimposed 75 3.5 Predicted signal and pointwise confidence bands joined by straight lines of incidence of melanoma data using an L-spline and a polynomial spline of order 7 80 3.6 Predicted innovation residuals of incidence of melanoma data with error bands for the L-spline and polynomial spline of order 7 models 81 3.7 Predicted signal and pointwise confidence bands joined by straight lines of incidence of melanoma data using a cubic spline and autocorrelation function of the innovation residuals 82 ix Acknowledgements I am most grateful to my research supervisor, Dr . Piet de Jong, whose advice and words were always a source of inspiration and guidance. I thank him for his patience and care to transmit not only the technical aspects of the statistical science but also the epistemological ones. Our extended discussions were most formative and have definitely marked these first steps of my carrer. I wish to thank Dr. Nancy Heckman, Dr. John Petkau, Dr. Ruben Zamar, and Dr. Harry Joe for the invaluable advice and help provided throughout my graduate studies. Thanks to Laurie Mac Will iam, Cheryl Garden, Matias Salibian, Birgitte R0nn, and Felix Menden who, as friends and/or collegues provided at different stages of my graduate student life uninterested support, love and care. I would also like to thank Piet de Jong for his financial support through an N S E R C grant. Finally, I want to thank my family. Without their love and trust this would have never been accomplished. Chapter 1 Introduction 1.1 M o t i v a t i o n Suppose that pairs of data values (yi,Xi), i = 1 , . . . , n , are observed and one wishes to describe the relationship between the components of the pairs. This is called a s m o o t h i n g p r o b l e m . The reason for the name is that one would like to describe the relationship between the variables without taking into account too much of the inherent noise which is present in almost all data sets. For example, if the x's are the times (or dates) at which the ?/'s were recorded, the set of pairs ( ? / ; , £ ; ) , i = l , . . . , n , is denoted a time series and the interest lies in describing how y evolves in time. A great deal of effort has been put by researchers to study time series and statistical models and inference tools have been devised for this purpose, for example, A R I M A models (see Box, Jenkins and Reinsel, 1994), the more modern scope of State Space models (see for example Anderson and Moore, 1979, Harvey, 1989 and 1993, West and Harrison, 1989) and in particular Structural T i m e Series models (see Harvey, 1989 and 1993). A l l these studies on time series can be translated with few modifications to the case where x is not necessarily time but any other equally spaced ordering variable for the ?/'s. In this thesis no assumption 1 on the spacing of the indexing variable, x, will be made. Another area of statistics where a lot of effort has been put into studying smoothing problems is nonparametric regression. In this case, the relationship between the variables is assumed to have an unspecified functional form g such that yi = g(xi) + ei, i = l , . . . , n , (1.1) where the e;'s are i . i . d . error terms. The main idea behind nonparametric regression is, in contrast to parametric models, not to make any assumptions on the form of g. In parametric models g would be known, at least implicitly, up to some parameters (e.g. polynomials). In non-parametric regression, qualitative assumptions on g, commonly the continuity of some derivative, are made in order to restrict the search of a suitable regression function to a reasonably small set. There are different approaches to estimate g in a non-parametric context. Perhaps the most popular non-parametric approach is the one based on polynomial splines. The estimation criterion that gives rise to polynomial splines consists on finding a function g that has certain regularity conditions (e.g. a degree of differentiability, square integrability of the derivatives, etc.), fits the data well and has some degree of smoothness (where smoothness or curvature is measured as the integral of the square of some derivative of g). It turns out that the polynomial spline of order 2m — 1 is a smooth piecewise polynomial of order 2m — 1 (see Chapter 4). It is claimed that because of the segmented nature of piecewise polynomials, polynomial smoothing splines will estimate the regression function with high fidelity and smoothness. It has been realised, though, that polynomial smoothing splines do not always provide satisfactory regression estimates. As a way of generalizing smoothing splines and at the same time cover the cases where smoothing polynomial splines fail to perform appropriately, attention is being paid to the so called 2 L-splines (see Chapter 3 for a definition). Most of the work done in the area of spline smoothing is not based on the use of stochastic models, but on the use of functional analysis theory. It is the claim of statisti-cians advocating nonparametric regression that by only supplying qualitative information about the regression function the data will speak for themselves concerning the actual form of the regression curve. As will be shown, parametric and nonparametric techniques do not cover the range of all possible approaches to solving smoothing problems. From the literature in non-parametric regression (e.g. Eubank, 1988) it seems as though the dichotomy "parametric vs. nonparametric" techniques for regression comprise the universe of possible techniques so far developed. There are other approaches, based on stochastic modeling, which have coexisted most of the time with the deterministic modeling techniques but have been somewhat relegated to equally spaced time series analysis and ignored by the researchers working in other regression (or smoothing) problems. The most important objective in data analysis is to understand the true law that dictates the behaviour of observed variables. For this purpose, usually the researcher is presented with not much more than observed data and little or almost nonexistent theory inherent to the variables for this purpose. In this case a sensible paradigm to adopt would be as follows. i) Recognize important characteristics of the data. ii) Develop a way to analyse the data, i.e. a methodology for data analysis and infer-ence, so that this adopted methodology can incorporate as much available perception 3 and theoretical background as possible. iii) After the data have been analysed according to step ii), check that the adopted methodology has captured the features in i) (by means of diagnostic checks) and look for further improvement if necessary. One can represent the above paradigm schematically as in Figure 1.1. i) recognise characteristics of data ii) adopt methodology for analysis data analysis iii) verify methodology is appropriate diagnostic checks Figure 1.1: A good data analysis paradigm. It is important that the choices of methodologies in step ii) are abundant so that the features of the data recognised in step i) can be captured. These methodologies have to be not only abundant but also easy to implement in terms of the feasibility, amount and level of complexity of computation involved. Also the methodologies of ii) should allow step iii) of the paradigm so that the loop can continue until a satisfactory description of the data is obtained. This type of paradigm, or the chart depicted in Figure 1.1 is not new in the literature. Box has presented similar schemes in his works (see for example Box, 1979, 1980, and Box, Jenkins and Reinsel, 1994 p.17). What is different is that Box does not make explicit step i) of the model building process and starts up in step ii). Other authors, inspired by Box's model building charts, have also designed and proposed more complete schemes 4 (see Cook and Weisberg, 1982, p.7, W i l d , 1994, Tong, 1990). For many regression problems it has been shown in the existing literature that the above paradigm is realizable by adopting in step ii) stochastic models such as state space models (see Chapter 2 for a definition and properties). The adoption of state space models permits the incorporation of observed features of data sets (for example in the form of components, like trends, seasonal and/or cyclic components, etc.) as well as existing background theory on the relationship of the variables (such as influence of covariates and the fulfillment of differential equations) in a specific and transparent way. Estimation within the context of state space models is relatively easy and convenient algorithms such as the Kalman filter, diffuse Kalman filter and smoothing algorithms are nowadays easily implemented on practically any computer. Diagnostic tools and tests are readily available to perform part iii) of the paradigm. O n the other hand, if one wishes to adopt the above data analysis paradigm and one decides to use polynomial spline smoothing as a methodology for data analysis, two problems arise. The first one is that in no obvious way is it clear how to incorporate observed features or underlying theory on the behaviour of the data in the data analysis methodology. The second problem is that because of the lack of diagnostic tools and the lack of a model the completion of the cycle indicated in Figure 1.1 cannot be done. As will be shown polynomial spline smoothing is a rather narrow class of methodologies to adopt for data analysis. The amount and level of complexity of computations involved in spline methods is not to be overlooked. Only for the case of cubic polynomial splines do there exist relatively easy algorithms and methods of computation (see Green and Silverman, 1994). Examples of articles where it is shown that smoothing by means of splines (either 5 polynomial or L-splines) is a problem that can be posed as a signal extraction problem using a particular stochastic representation of the unobserved signal or, in other words, as a smoothing problem in the context of a particular state space model are those by Weinert, Byrd and Sidhu (1980), Wecker and Ansley (1983), Ansley and K o h n (1986), Kohn and Ansley (1987), and Ansley, Kohn and Wong (1993). Weinert, Byrd and Sidhu (1980) studied the relationship between L-splines and stochastic processes and gave a state space model to compute L-splines. Wecker and Ansley (1983), found a state-space model for polynomial splines. Ansley and Kohn (1986) and K o h n and Ansley (1987), and Ansley, Kohn and Wong (1993), among others, based on Weinert, Byrd and Sidhu's work of 1980, explore further properties and extensions of the relationship between L-splines and stochastic processes. As Box (1979) points out, it should be recognised that the choice of a methodology for data analysis is not an irrevocable decision and the investigator need not attempt to allow for all contingencies a priori. A good data analysis methodology should allow for what Box calls a criticism stage that allows understanding of an initial analysis. This should lead the researcher to question the methodology adopted or to look for improvement of the initial analysis. For all the reasons above, it is proposed in this thesis to adopt the paradigm repre-sented in Figure 1.1 for regression analysis. 1.2 C o n t r i b u t i o n s In this thesis, the problem of smoothing using nonparametric tools and stochastic models is addressed. The nonparametric tools considered in this thesis are polynomial smoothing 6 splines and L-splines and the stochastic models are state space models. The relationship of both approaches, nonparametric and stochastic models, is studied. A comparison of the performance of these methods is done by analysing real data. In order to carry out the comparison of these two approaches, continuous state space models are studied in detail. First, it is studied how to derive a discrete state space model from a continuous state model. The most difficult part in this derivation is to compute the state error covariance matrix. When the stochastic process defined by the state equation of the continuous model is stationary then a result exists that makes the computation of the matrix easy. The problem when the process is nonstationary has so far not been addressed in the literature. This problem has importance since most of the stochastic processes that are associated with smoothing splines (L- and polynomial) are nonstationary. A result, applicable to the nonstationary case, is shown that makes the computation of the state error covariance easier than the direct calculation, in many practical cases. Also, it is shown how to initialise the diffuse Kalman filter to carry out estimation and prediction in the context of a discrete model derived from the continuous one. In general it is shown how to translate all the existing estimation and prediction tools for the discrete case to the continuous case. Smoothing by minimising a penalised least squares criterion is associated with smooth-ing by using a certain continuous state space model. This state space model, in many cases, is not always properly specified in that the initial conditions are always diffuse. B y means of the analysis of an example it is shown that proper initial conditions should not be overlooked since they affect in great extent the estimation and prediction results. It is also shown, in this same example and a second one, the advantages of a model based approach to smoothing. It is shown how different models can be compared and how the 7 paradigm in Figure 1.1 can be realised using state space models. In order to compare polynomial spline smoothing with a stochastic model based ap-proach to smoothing, the penalised least squares criterion giving rise to polynomial splines is studied in detail. In particular it is shown that the smoothing polynomial spline defining criterion is equivalent to another penalised least squares criterion. This result generalises a known result for cubic splines. This new penalised least squares criterion allows to make the properties of smoothing polynomial splines transparent. In particular it is shown that a stochastic model associated with smoothing polynomial splines is an A R I M A plus noise model where the only parameter to be estimated from the data is the noise variance. 1.3 O u t l i n e The rest of this thesis is organized as follows. Chapter 2 is a review of state space models, where the main ideas on state space modeling, definitions, estimation tools and results are presented. This chapter should bring ideas of how these stochastic local models can be conceived and how to carry out estimation and prediction in the state space model context. Also, this chapter is fundamental for the development of the next Chapter. In Chapter 3, estimation within the continuous state space model is addressed. New results on the implementation of estimation and prediction tools (diffuse Kalman filter and diffuse smoothing) are presented. In particular, initialisation issues of the diffuse Kalman filter are studied. Also, L-splines, a nonparametric smoothing tool, are briefly reviewed. The focus is to view L-splines in the context of continuous state space models. Examples are given to illustrate how nonparametric regression methods can be easily implemented 8 when using a model based approach. The advantages of the model based approach to nonparametric smoothing are illustrated with examples. It is shown how the exploration of a model, rather than the construction and minimisation of an objective function, can result in a fruitful endeavour. In Chapter 4 , a review of polynomial smoothing splines is given. Mainly, a basis for the space of smoothing polynomial splines (of any pre-fixed order) is studied so as to extract and make more transparent certain properties of polynomial smoothing splines. It is shown that polynomial smoothing splines can be obtained by minimising a certain penalised least squares criterion, generalising a well known result so far proved for cubic splines only. In Chapter 5 polynomial splines are further explored and a connection between poly-nomial splines, solutions to penalised least squares problems and A R I M A processes is shown. In particular it is shown that polynomial smoothing of order 2m — 1 is equivalent to A R I M A ( 0 , m, m — 1) plus noise modeling where all the M A parameters are fixed. This result may bring some more light in the understanding of the polynomial spline smoothing tool. 9 Chapte r 2 The State space mode l A general class of models of much current interest is that of state space models, also referred to as S S M from now on. State space models were originally developed by control engineers, particularly for applications concerning navigation systems such as controlling the position of a space rocket. Currently, and especially due to the contributions of Anderson and Moore (1979), Harvey (e.g. 1989), West and Harrison (1989), de Jong (e.g. 1987,1988,1991), among others, S S M have been found useful in time series problems. It is the purpose of this thesis to show that S S M can be useful in general regression problems as well. Suppose that pairs of data values (yi,X{), i — 1 , . . . ,n, are observed. The smoothing or regression problem consists of describing the relationship between the components of the pairs. B y this notation, the variable xt- plays the role of an indexing v a r i a b l e . That is the y's are ordered according to the x's and x\ < • • • ,xn. If x is time, then the data are denoted a time series and the interest lies in describing the evolution of y in time. If y is energy consumption and x is temperature then the interest lies in the study of the changes of energy consumption as the temperature changes. If y is the height of 10 a child and x its age in months then the interest is to study how the child grows in height as it grows older. The indexing variable plays the role of imparting structure to the data; it defines contiguity of the y's. Think, for example, of the y's being energy consumption of a household. Ordering the y's in the (time) sequence they were observed usually gives a completely different structure to the data than if the j/'s had been ordered by a temperature indexing variable. When the researcher measures any sort of (unobservable) signal, it will typically be disturbed by noise so that the actual observation is given by observation=signal-(-noise. In state space models the signal is taken to be a linear combination of a set of variables, which constitute what is called the state vector. The state vector at X{ denotes the so called "state of nature" at x ± . Note that the state vector is unobservable and that it has to be predicted from the data. Often the state has a physical meaning. For example, the progress of a spacecraft can be monitored using a S S M with a state whose components describe the velocity, acceleration, coordinates, and rate of fuel consumption of the spacecraft. The N A S A Space Program utilises S S M to control its spacecrafts. In many economic contexts the state may represent constructs such as trends or seasonalities. For example, Crafts et.al. (1989) employ the S S M to predict the historical unobserved trend and cycle components of the British industrial production index. In other situations, where physical information is not available, the state can represent an unobserved and/or unknown process driving the observed system. For example Watson and Engle (1983) use the S S M to estimate unobserved wage rates. Other applications of the S S M include a study in inventory 11 control (Downing et.al., 1980) and of statistical quality control (Phadke, 1981). In the area of policy, Harvey and Durbin (1986) conducted an interesting study of the impact of seat belt legislation on road casualties. In parametric regression the state vector stays constant as the indexing variable, x, varies and the constants that enter in the linear combination of state variables change. In both, parametric and nonparametric regression, the signal is assumed to be de-terministic. In parametric regression the functional form of the signal is assumed to be completely known up to some parameters. In nonparametric regression the signal is not assumed to have a specific form but to have some general qualities, such as a certain degree of differentiability. It is usual to read, especially in the nonparametric literature, that parametric vs. non-parametric regression constitute an ambivalence that covers in great extension a universe of alternatives for the solution of regression problems. As will be seen throughout this chapter and the following ones, there exist other options to solve regression problems. A great deal of effort has been put into studying signal plus noise models where the signal is stochastic, that is the signal is assumed to be a stochastic process rather than a deter-ministic function. Hence the way that the signal at Xi+\ is related to the signal at Xi is established in a probabilistic rather than in a functional manner. The study of the signal as a stochastic process has yielded interesting results that allow insightful solutions of the problem of estimating the unknown signal. Perhaps because these studies were mo-tivated by econometric applications, they were relegated to the study of equally spaced time series. Typically, and this is a strength of state space models, the state vector changes as the 12 indexing variable, x, changes (whereas the constants that enter in the linear combination of the components of the state vector to give the signal are usually fixed.) This is because the state is modeled as a stochastic process on x that suits, or is thought to suit the data well. Hence, S S M are stochastic local models for the data where local refers to the adjacency of the data dictated by the ordering variable, x. As a consequence, S S M are robust in the sense that isolated departures from the bulk of data will have a limited influence in the global results of estimation. In this way, smooth estimates of the signal that produce a good local representation of the data can be obtained. As an example, suppose that the ordering variable, x, is equally spaced, with a;,-_i—x; = 1, and that it is believed that the relationship between y and x is linear. Then, yi = a + b x i , i — 1 , . . . , n, for some a, b G R. This is equivalent to saying that = yi + b, with starting condition yo = a. If one wants a more flexible model than the initial linear model for the relationship between y and x, then, starting from the linear model one could add a disturbance, e, and postulate that yi = a + bxi + e,-. Or, one could add a disturbance, 77, to the alternative formulation and say that ?/;+,• = y,- + b + rji. These two new "more flexible" formulations are very different. The first one, yi — a + bxi + et'> is the usual linear regression paradigm, whose drawbacks, for example lack of flexibility or lack of robustness, are well known. The second one, j/,-+t- = yi + 6 + rji, is a state space description that produces a l o c a l l i n e a r d e s c r i p t i o n of the data, it represents local characteristics of the data. Therefore, an initial statement of a model and the placement of disturbances is crucial in making this initial model more flexible. The construction of S S M is concerned in identifying possible components (the state) of an initial global formulation of the signal and the placement of random disturbances in the components in order to obtain a local description. 13 Following, a review of discrete state space models is presented. Discrete state space models have been thoroughly studied especially to use these models in econometric ap-plications, e.g. time series. As will be seen, these models can be used in most regression situations as well. The formulation and treatment of discrete S S M involves keeping track of the particular spacing of the indexing or ordering variable, x. The next chapter intro-duces continuous state space models which are formulated independently of any particular spacing of a;'s. Once the data are collected, a discrete model is derived from the contin-uous one in order to suit the data. Most of the results derived for the discrete case will apply fairly directly. In light of the above paragraph, the review of results on discrete S S M that follow is essential to the development of Chapter 3. In section 2.1 a formal definition of S S M is given. Section 2.2 shows some examples and explores in them the usefulness of the formulation of S S M in certain situations. The rest of the chapter deals with technical details about the implementation of the prediction and estimation tools used with the S S M . Section 2.3 deals with the criteria that are used for estimation in the context of S S M and the Kalman and smoothing filters which are the tools to be used to implement the estimation. Sections 2.4 and 2.5 deal with diffuse S S M , and S S M which have applied since time immemorial. These models are useful generalisations of the original S S M of Section 2.1. Section 2.6 shows how generalised cross validation ideas can be implemented in the context of S S M if a distribution free method for estimation is preferred. 2.1 D e f i n i t i o n o f the state space m o d e l State space models offer the possibility of constructing useful local models to solve smooth-ing problems. The advantage of explicit stochastic local models is that they are more 14 flexible than parametric regression models. In S S M the data play an important role in determining the local changes. This brings also the advantage of having robust models to local departures from the bulk of the data. Also, important observed features of the data, such as trends or cycles, may be incorporated in the model. In this section the following notation is used: • If a is a random vector then Var(a) denotes its variance-covariance matrix. • 7 ~ (c, a2C) indicates that 7 is a random vector with £ ' ( 7 ) = c and Var(-y) = <r2C. Expectations and covariances are always unconditional. def • If A and B have identical column dimension, then (A; B) = (A1, B')'. • r/* indicates the number of components in y. The following formulation of the state space model is the one given by de Jong (1991, and thereafter). It is somewhat different from the one presented in Harvey (1989) or West and Harrison (1989). This formulation is more general than the other existing formulations. It also allows, in an elegant way, an extension of this model to the diffuse state space model. D e f i n i t i o n 2.1 A random vector y = (y\, y2;...; yn) is said to be generated by a state space model, denoted y <^SSM, if yi - Xi/3 + ZiCti + GiUi, i = 1,..., n; (2.1) Qfi+i — Wif3 + T ;CY; + H{Ui, i = 0 . . . ,71 (2.2) where 15 (0) Xi, Zi, Gi, Wi, Ti, and Hi art known matrices. (1) (u0; U\\...; un) ~ (0, a2I) with a2 > 0, (ii) cto = 0 and = b + Bj, where 7 ~ (c, o~2C), b is fixed and B has full column rank, (iii) 7 and (UQ; ui;...; un) are uncorrelated, (iv) C is non-singular unless C = 0; Var(y) is non-singular if C = 0. It is always the case that pairs (yi, x t), i — 1 , . . . , n, are observed, where the x's play the role of an ordering or indexing variable, hence x\ < . . . < xn. If the z's are the times at which the y's were observed then we are dealing with time series data. T i m e series are usually observed at equally spaced intervals. Associated with the order of the ?/'s is the notion of proximity of the ?/'s. We say that yi and yj are contiguous or close if Xi is close in Euclidean distance to Xj. This notion of closeness of the y's is important since we are dealing with local models. Hence, in time series, closeness of the ?/'s means that the y's were observed relatively close in time. Note that the y's can be multidimensional. Equation 2.1 is called the measurement or observation equation. It shows how the measurements or observations are generated. The signal at Xi is Xifi + ZT,Q;,- and the measurement error at Xi is GiUi. The vector a is called the state vector and on = a(xi) is called the state at Xi. Equation 2.2 is the state or transition equation. It shows the particular stochastic process, or statistical dynamics, of the state, a. The quantity HiUi will be referred to as the state error. Often, if the equally spaced then Ti = T, that is, the state transition matrix, Ti, is constant and does not change as x changes. This means that the dynamics of the state are constant over each equidistant interval of the a?'s. 16 If C = 0, note that there is no loss of generality in assuming that cto = 0. If «o = 0 then E(ai) = W0/3 and Var(a-C) = O-2H0HQ. If in fact the data being modeled dictate that E(a0) = a and Var(a0) = a2AA', where a 0 is uncorrelated with (u 0 ; i t i ; . . . ; u n ) then e*i = (Wo, r0)(/3; a) + {H0, A)(u0; u_ i ) , where £ ( u i ) = 0 and V a r ( u i ) = a2I. Redefining f3 as (/?;a) shows that the above formulation is equivalent to assuming Q.Q = 0. Neither is there loss of generality in assuming that f3 and U{ are the same in the observation and state equations. In particular through appropriate choice of Xi, Wi, Gi and Hi distinct components of /3 and u,- can be included in or excluded from either equation. Different specifications of the S S M of special interest are: • C — 0 and so (3 is fixed; • C is non-singular and /3 is random; • ~ * 0 a n d P is said to be diffuse. These different specifications of the S S M will be dealt with separately in the following sections. 17 2.2 S o m e e x a m p l e s of state space m o d e l s 2.2.1 The parametric regression model as a S S M Consider the case where pairs (?/,-,x;), i = 1 , . . . , n , are observed and it is assumed that yi = Si + ti where Si is the signal at Xi and e; is a (0,cr2) random variable. Further it is assumed that s,- = g(xi,f3) where g is a known function and (3 is generally an unknown vector of parameters. In the case of linear regression, g(x{,/3) = X(xi)f3, where X(x,) is observed together with the pairs (yi, X j ) , i = 1 , . . . , n. This is a trivial kind of S S M with Xi = X(xi), Zi = 0, Gi = 1. The state equation, and hence the starting condition ao, is redundant. Hence, there are no dynamics in the linear regression model. The specification f3 — b + Bj allows for both fixed and random effects in f3. If (3 is known then B = 0. Heteroskedasticity can be dealt with by appropriate definition of Gi. Another possible way of stating the linear regression model in state space form is to set Xi = 0, Zi = X(xi), G,- = 1, Wi = 0, = / , Hi = 0 and hence the state is deterministic and constant. In the case of polynomial regression, g(xi, (3) = (1, Xi, x2,..., x™ )/3 and Xi = (1, X{, x2,..., x™), Zi = 0, Gi = 1. Again, the state equation is in this case superfluous. The S S M contains the linear and polynomial regression models as extreme cases, where there are no dynamics. 18 2.2.2 Structural models In general, and in particular in this thesis, this is how S S M will be constructed: taking a particular deterministic model for the signal, studying its linear dynamic properties and adding stochastic disturbances to the deterministic model appropriately so that a local statistical description is obtained. To illustrate how this is performed in practice, the example of structural models will be reviewed. Structural models are a special case of state space models. These models were fully studied by Harvey (e.g. 1989), especially in the equally spaced time series context. Har-vey and Stock (1986) and Harvey (1989, C h . 9) also study structural time series models in continuous time. This approach, as already shown in Harvey and Stock (1986), and as reviewed in Chapter 3, provides a way to handle the cases of unequally spaced observa-tions. In this section only structural models stated in discrete time will be reviewed. In Chapter 3 the continuous formulation of some structural models will be given. No specific assumption about the spacing of the ordering variable x will be made. A structural model consists of a signal plus noise model for the observed variable y, indexed by the variable x, namely yi = Si + et, i = 1 , . . . ,n (2.3) where yi is the observation of y made at Xi, 6t' is a (0, cr2) random variable independent of 5,-, and Si is an unobserved signal. In a structural model, the signal can be written as a sum of different components which have a direct interpretation in terms of the observed data. Also, these components are assumed to have a specific form, usually termed "linear dynamic". 19 As an example, one may write Si — fli -{~ 7t where \ii is the so called trend and -ji is a seasonal component. To illustrate the meaning of a trend and seasonal component one may think of the sales of a department store observed at different periods of time. It is known that sales increase during certain times level of sales of the store. Each of these components has a significance and impact of its own and it is important to be able to separate them in order to draw conclusions from the data. As already said, it is assumed that each of the components has a "linear dynamic" representation. W i t h this we mean that a component at Xi is a linear function of the component at plus some disturbance. In the case of the trend, fi, suppose that we have reason to believe that it is approximately a linear function of x. Hence, if fi were perfectly linear one would have that fj,(x) = a + bx for some a and b. Therefore, fi(xi) — fi(xi-i) = bhi, with hi = — X{ and so, taking /i,- = p{xi) and /?,• = b, one can write Deterministic models for the trend, like (2.4), would in general produce a poor de-scription of "noisy" data. Hence a generalization, for example, of the year (e.g. Christmas) and this introduces a seasonal effect in the actual trend or (2.4) f i i + 1 - m + fSihi + rji Pi+l = Pi + Cii (2.5) where rji and Q denote (0, a2.) and (0, cr2.) disturbance variables respectively, can result in a more appropriate or flexible description of the trend. Here, cr2. and Q2 depend on hi since 20 closer y's will have more in common than further apart ones and hence the disturbance in a narrow interval will be smaller. If the IT'S are equally spaced then the 77;'s and £ ' s have constant variances. If CTQ > 0 then we allow the slope of the trend to vary in each of the intervals determined by the x's. So, one can appreciate how, starting from an idealised and global deterministic signal and by adding stochastic disturbances appropriately one can construct a flexible stochas-tic model which will focus on local characteristics of the data. Adding disturbances to components is more effective than just adding a single disturbance to the deterministic signal. Note how a completely deterministic treatment of the signal, as in the previous section, in the case of parametric regression, does not suggest other possible ways of mak-ing the model more flexible except for turning a parametric signal into a more unspecified one. For the seasonal component, one may assume that it is approximately a trigonometric function of x. If it were perfectly trigonometric then it would be of the form j(x) = asin(cux) + bcos(u>x) for some a , b, and u>, where w is a frequency corresponding to period p , that is to = 2 i r / p . This ensures 7(3;) = 7 ( 2 + p). Then, defining -f*(x) = a c o s ( u x ) — bsm(u>x), one has that j ( x i + 1 ) = 7(2,- + hi) = a s i n ( u j ( x i + h i ) ) + b c o s ( u ( x i + h i ) ) — [a s'm(ujXi) + 6cos(a;a;j)] c o s ( u h i ) + [a cos(ujXi) — 6sin(a;a;j-)] sm(u>hi) and 7 * ( x i + i ) = y * ( x i + h i ) = a c o s ( u ( x i + h i ) ) — bsm(uj(xi + h i ) ) 21 = [acos(ujXi) — 6sin(u;a;,-)] cos(u>hi) — [ a s i n ( a ; X j ) + bcos(uxi)] sm(uhi). Letting 7,- = 7(2:;) and 7* = j * ( x i ) one has that ( 7*'+i \ — c o s ( u h i ) s m ( u h i ) \ 7 * + i / _ - s m ( u h i ) c o s ( u h i ) ( ) (2-6) Since in real life one is unlikely to observe an exactly trigonometric seasonal signal, one may add disturbances in (2.6) so that depend on h{. In this way, once one has a dynamic linear representation for an idealised signal, one may construct a more flexible model by adding disturbances appropriately. Following are some of the most common structural models used in practice. P u r e l o c a l l inear t r e n d m o d e l In this case we have a model where only a trend is observed. This trend is assumed to be locally linear with added disturbances. The model expressed in state space form would be where the e;'s are independent of the 77^ 's and £-'s e^ 's are i . i .d . (0,cr2), and 7/t- and are Table 2.1 shows how the placement or not of disturbances gives different model spec-ifications for the trend. I 7*'+i j _ COS(UJhi) s i n ( u h i ) \ 7 * + i / _ s m ( u h i ) cos(wftj) where K ; and K* are (0, A;,) and (0, k*) variables, re (2.8) (0 , (J 2 . ) a n d (0,<7 2.) r e s p e c t i v e l y w i t h c o v a r i a n c e 22 Type of trend model < ° l locally linear > 0 > 0 stochastic (or local) level with fixed slope > 0 = 0 smooth trend = 0 > 0 deterministic linear = 0 = 0 Table 2.1: Different trend models. Autoregressive component plus noise model In (2.5) let //,• = 0 for all i and, moreover, assume that = pf3{-\ + 0 < p < 1. The signal then follows an autoregressive process of order 1. If p = 1, then the signal follows a random walk process. The corresponding state space formulation of the model is Vi — OL{ -\- t{ a i + 1 = poti + Ci where e,-'s and the £,-'s are independent, e;'s are i . i .d (0, a2), and the £,'s are (0, a2.), where a2, depends on h{. This model is frequently used in absence of physical information about the signal. Trigonometric cycle or seasonality model Since seasonality is a periodic phenomenon it is natural to model it using trigonometric functions. 23 The corresponding state space formulation is yi = (l_0)a t- + e,-c o s ( u h i ) s m ( u h i ) s'm(uhi) c o s ( u h i ) a,+ ;: , ™ where K i and K* are correlated white noise disturbances, usually assumed to have a com-mon variance a2.., with a2., depending on hi, and p is a "damping factor" which lies in the range 0 < p < 1. There are further variations of these models and only some of them were illustrated. 2.2.3 A R I M A models as state space models A complete reference for A R I M A models can be found in Box, Jenkins and Reinsel (1994). A R I M A models are based on the idea that a time series, yi, i = 1 ,2 ,3 , . . . , in which successive values are highly dependent, can frequently be regarded as generated from a series of independent (0,a2) "shocks" or disturbances The process ut-, u ; -2 , • • • is called a discrete white noise process. It is assumed that some sort of "linear filter" transforms the white noise process into the stochastic process of interest yielding yi = p, + Ui + ipiUi-i + f/>2^t-2 + • • where p represents the level of the process. The attention is focused on modeling the signal as a stationary process, an assumption which is crucial for estimation and modeling purposes. The above idea is conceptually based on a result by Wold (1938), who established that any zero-mean purely non-deterministic stationary process ?/,• possesses a linear representa-tion yi = Ui + YsJLi V'j^i-j) with the u,'s being (0, a2) and uncorrelated, and Y^JLi ip] < °°-In this sense, the models that are studied are moving average ( M A ) models, where the process j/t- is assumed to be a finite linear combination of i . i .d . disturbances, and 24 the autoregressive models (AR) , where yi is assumed to be a finite linear combination of past values of yi plus a disturbance. Back substitution in the past values of j/j allows expressing an A R process as an infinite linear combination of iid disturbances. To achieve greater flexibility in fitting an actual time series, yi can also be assumed to be the sum of an autoregressive process and a moving average process, or in the Box and Jenkins terminology, an A R M A process. Box and Jenkins recognise the importance of having at hand some nonstationary models, in addition to stationary A R M A models, since some real life processes in industry, business, etc., are better represented as such. In particular they consider a special type of nonstationarity. They assume that if the process is nonstationary then, after differencing it a few times it will have a stationary A R M A representation. These type of nonstationary models are called integrated moving average models or A R I M A , for short. In using A R I M A models, it is assumed that the values of the indexing variable, x, are equally spaced. Also, let us define the backward shift operator, B, by B3yi = yi-j, This operator is extensively used in the A R I M A literature. A formal definition of the mentioned processes follows. Def in i t ion 2.2 The process yi is said to follow an ARMA(p,q) model if yi = 8 + ai?/i_i + a2yi-2 H h apy^p + et- + &ie,-_i H h &,e t-_„ where e,- ~ (0, cr2I) and yo, y~\,..., y~v are unobserved and may be either fixed or random. Also, av ^ 0, and bq ^ 0. The process yi is said to follow an ARIMA(p,d,q) model if (1 — B)dyi follows an ARMA(p,q) model. 25 As a consequence of the definition, if yi follows an A R I M A ( p , d, q) model it also follows an A R M A ( m a x ( p , d), q) (not necessarily stationary) model . It turns out that A R M A models can be cast into S S M form, as the following result shows. T h e o r e m 2.1 (Gardner et.al., 1980.) Suppose that yi follows an ARMA(p,q) model. Put c = max(p, q + 1) and define T = { (a i ; . . . ; ac), (7; 0)}, H = (1; b\;...; 6 c_i) and W = (1; 0 ; . . . ; 0) where the aj or bj for respectively j > p and j > q are defined as zero. Then, yi — ( l ,0, . . . ,0)a,-a,-+i = Wp3 + Tai + Hui where u. (0 ,a 2 ) . E x a m p l e 2.1 (ARIMA(0,2,1) models as SSM.) Let yi be generated by an A R I M A ( 0 , 2 , 1 ) process. Then yi - 2j/,-_i + yi-2 = Ui + where u8's are independent and it; (0 ,<7 2 / ) . Then, yi can be generated from the following state space model yi = (1 0)ai 26 It is easy to see that if P = 1 0 -1 0 then P-^TP = 1 1 0 1 Hence the ARIMA(0,2 ,1) S S M is equivalent to V, = (1 OK, Note the similarity with the local linear trend model (2.8). The state transition ma-trix is the same, but there is a different placement of disturbances. In particular, the measurement error is not present. In this way one can view a local linear trend model as an "ARIMA(0,2 ,1) plus noise" model. • Some advantages of casting A R I M A models as S S M are the following. • State space models permit the treatment of A R I M A models in a unified framework. • The A R M A modeling and estimation theory relies heavily on assuming stationary stochastic models for the signal, or at least models that become stationary after taking a finite number of differences. Perhaps the stationarity concept is better suited for very long time series where an estimated "equilibrium" probability dis-tribution could be justified. In many regression problems it is not always the case that the indexing variable, x, can be expected to take values in arbitrarily large intervals, unlike most time series situations. Think for example of studying energy consumption as a function of a climatological variable such as temperature. It is unlikely that temperatures exceed a certain finite range. O n the other hand raw 27 nonstationary components may be of more interest than a transformation of the component to achieve stationarity. • The equal spacing of the indexing variable, x, is a big restriction. Many regression problems have an unequal spacing of the variable. B y assuming many missing observations, one could still try fitting A R M A models but the implementation of the estimation procedure would become difficult and cumbersome. As shown above, A R I M A models can be seen as S S M and hence S S M are more general. As it has already been shown, the nonstationarity issues (e.g. de Jong and Chu-Chun-L i n , 1992) and unequal spacing (e.g. Harvey and Stock, 1986) can be circumvented in some cases in the S S M context. In this thesis more issues on nonstationarity and unequal spacing in the S S M will be dealt with. 2.3 P r i n c i p l e s a n d tools o f p r e d i c t i o n a n d e s t i m a -t i o n for the S S M In using observed data to infer about random variables, the most common problems are the following. • Forecasting: given ( y i , x - i ) , . . . , ( y „ , x n ) it is desired to obtain a prediction of y n + k at x n + k , k > 0. • Interpolation: given (j/ l 5 x x ) , . . . , (?/;_!, x ^ ) , ( y i + 1 , x i + i ) , . . . , ( y n , x n ) one wants to predict yt- at X{. • Smoothing: if the observations are generated by a S S M , a prediction of the unob-served state a{ at Xi is desired. Another smoothing problem is the signal extraction 28 problem where a prediction of the signal Xid + zT,o:,- at X{ is desired. Also, an estimate of the measurement and/or state error components are sometimes desired. Frequently in the S S M there are unknown parameters, for example some disturbance variances and covariances. Prediction and smoothing can only be carried out once these parameters have been estimated. Some ways to proceed with the estimation of these parameters is via maximum likelihood or via cross-validation. It will be shown that the same tool that will be used for prediction can be used to obtain an expression for the likelihood of the j/'s and/or an expression for the cross-validation criterion. Maximum likelihood is to be used when the disturbances, u,'s, are assumed to be normally dis-tributed. If the normality of the disturbances is doubted then a distribution-free method like cross-validation can be used. The prediction criterion to be used is the Mean square error (Mse) and the estimators that will be considered are a constant plus a linear combination of the observations. Henceforth, predictors will denote "best linear predictors" where "best" is relative to minimization of the Mse. Best linear prediction and minimum mean squared error are the key words in the implementation of the prediction procedures in the context of the S S M . The tools to be used are the Kalman and smoothing filters. These algorithms permit the incorporation of information sequentially which makes prediction and estimation methods easy to im-plement. In prediction procedures where all the data are incorporated at once, as in classical linear regression, quite cumbersome operations (for example the inversion of a large matrix) may have to be performed. 29 In this section definitions and technical aspects of inference in the context of the S S M , as per Definition 2.1, will be given. The proofs of the results will be omitted and can be found in textbooks like Anderson and Moore (1979). First, in Subsection 2.3.1 the concept and some properties of best linear prediction are formally stated. Subsections 2.3.2 and 2.3.3 contain some important properties of best linear predictors in terms of combining information. This is a preparation and the basic material for the Kalman filter algorithm, which is explained in Subsection 2.3.4. The Kalman filter provides predictions of the state at x;, given all the information collected at x i , . . . , x ; _ i . The smoothing filter gives predictions of the state given all the available information (fixed interval smoothing), a prediction of the state at a particular point xt-as the sample size increases (fixed point smoothing), and/or a prediction of the signal at a point Xi where no observation of y was made (interpolation). These smoothing situations are reviewed in Subsection 2.3.5. 2.3.1 Minimum mean square error prediction This subsection contains the definitions and some of the properties of best linear pre-dictors. This is the type of estimators that will be considered for estimation within the S S M . Suppose that X and Y are random vectors such that E(X) — \ix, E(Y) — \iy, Var(X) = a2Zx, Var(Y) = c r % , Cov(X,Y) = <J2Zxy. The prediction of X given Y will be the random variable denoted by P(X\Y). D e f i n i t i o n 2.3 The (best linear) predictor of X given Y is P(X\Y) = c + b'Y where c and b are such that E(X — c — b'Y)2 is minimum. Further, the mean square error of 30 prediction is Mse(X\Y) =J E (X - P(X\Y))2. The following theorem shows some important properties of the best linear predictor. Theorem 2.2 Suppose that X and Y are random vectors. Then, a) P(X\Y) = c + BY if and only if E(X - c - BY) = 0, Cov(X - BY, Y) = 0; b) ifTiy is nonsingular then P(X\Y) = fix + E ^ E ^ f y - py), Mse(X\Y) = a2 ( E X - E ^ E ^ E J ; c) IfY,xy = 0 then P(X\Y) = fix and Mse(X\Y) = a2Xx; d) Mse(X\Y) = Var(X - P(X\Y)) = Cov(X,X - P(X\Y)); e) P(X\Y) is unique. Theorem 2.2 asserts that the predictor is unbiased and the prediction error is always uncorrelated with the used variables. Also, if E Y is nonsingular then the predictor of X given Y is uniquely expressible. The condition Hxy = 0 is implied by E Y = 0 and hence the expressions in c) hold if E Y = 0. Part d) asserts that the mean square error matrix is in fact the variance matrix of the prediction error and gives simplified formulas for Mse(X\Y). It is worth noting that if (X; Y) is normally distributed then P(X\Y) — E(X\Y) and Mse(X\Y) = Var(X\Y), where E(X\Y) and Var(X\Y) are respectively the conditional mean and conditional variance of X given Y. 31 2.3.2 Combining predictors and iterated prediction The next result shows how predictors which use two uncorrelated pieces of information are combined. T h e o r e m 2.3 If X, Y and Z are random vectors, then P(X\Y,Z) = P(X\Z) + P(X\Y - P(Y\Z)) - E(X) Mse(X\Y,Z) = Mse(X\Z) + Mse(X\Y-P(Y\Z))-Var{X) T h e o r e m 2.4 If X, Y and Z are random vectors then P(X\Y) = P(P(X\Y,Z)\Y) Mse(X\Y) = Mse(X\Y,Z) + Mse(P(X\Y, Z)\Y) 2.3.3 The prediction and Kalman update formulas The results of this section relate to updating a prediction in the light of additional in-formation. This is important in situations where information arrives sequentially. In the case where all the information is available at once, then it can be incorporated piece by piece until all of it is used to get a prediction. In the theorems below, a is to be predicted and P(ct\z) is available, where Mse(a\z) = a2M. Additional information y = X(3+Za+Gu arrives, where /3 is fixed and u is random, with zero mean and uncorrelated with a and z and hence with P(a\z). Desired is P(a\y, z) and Mse(a\y, z). The usefulness of these results is in terms of the Kalman filter algorithm of the next section. 32 T h e o r e m 2.5 (Prediction update formula) Suppose that a and z are random vectors. Further suppose y = Xft + Zcx + Gu where ft is fixed, u ~ (0, a2I) and u is uncorrelated with a and z. Put P(a\z) = a, Mse(a\z) = a2M and B = MZ'(ZMZ' + GG')~X. Then, P(a\y, z) = a + B(y - Xft - Za), Mse(a\y, z) = a2(I - BZ)M. The next theorem deals with a slightly more general setup and is the key to the important algorithm known as Kalman Filter. T h e o r e m 2.6 (Kalman update formula). Suppose that a and z are random vectors. Further suppose y = Xft+Za-\-Gu where ft is fixed, u ~ (0,<r2/) andu is uncorrelated with a and z. Put P(a\z) = a, Mse{a\z) = a2M and K = (TMZ' + HG')(ZMZ' + GG')'1. Then, P(Wft + Ta + Hu\y,z) = Wft + Ta + K(y — Xft — Za), Mse(Wft + Ta + Hu\y,z) = a2[(T - KZ)MT' + (H - KG)H']. The previous theorem is basic to forecasting, signal extraction, smoothing and inter-polation. Usually, z is previous data, while y is the observation that is newly arrived. The matrix K is called the "Kalman gain" matrix. Note that actual predictions do not involve <r2 and that the matrices K and (T — KZ)MT' + (H — KG)H' do not involve z or y. Thus, these may be calculated before y is actually known. Note that a2D =7 a2(ZMZ' + GG') = Mse(y\z) and Mse(Wft + Ta + Hu\y,z) = a2[(T - KZ)MV + (H - KG)H'\ = a2(TMT + EE' - KDK') = Mse(Wft + Ta + Eu\z)-a2KDK'. 33 Hence a2KDK' is a measure of the extra information provided by y over and above that provided by z. 2.3.4 The Kalman filter The Kalman filter (KF) is an algorithm to recursively update the prediction of the state on in a S S M as additional observations come to hand. The following theorem uses the notation of Definition 2.1. Theorem 2.7 (Kalman Filter) Suppose that y <— SSM, with C = 0. Put a\ = WQ3, PL — HQH'Q. For i = 1, . . . , n, define e; = Yi — Xip3 — zTjOtj-, Di — ZiPiZ[ + GiG\, 1U = {TiPiZl + HiGfiDT1, ai+i - Wi/3 + Tiai + Kiei, Pi+1 = (Ti - KiZi)PiT[ + (Hi - KiGi)H'i. Then, ai+i = P{<xi+i\y\,---,yi) a2Pi+1 = M s e ( a i + l \ y 1 , y { ) e* = yi - P(yi\yi,...,yi-i) a2Di = Mse(yi\y1,...,yi-1). The quantities e,- are denoted "innovation errors" and if i ^ j, Cov(ei,ef) — 0. 34 The K F equations indicate how to recursively predict the state and the next observa-tion yi. However the K F has further ramifications for smoothing, k-step ahead prediction, interpolation and likelihood evaluation. Further points regarding the K F worth noting are: • The innovations are that part of yi not predictable from the previous values of the time series. If ?/,• is scalar then so are et- and Di. The theorem asserts that innovations are uncorrelated. • The matrix K is called the Kalman gain. Note how the "new" prediction of is a linear combination of the old prediction and the prediction error and there is no need to remember the previous data. • The K F is recursive. This is important in situations when observations arrive se-quentially in time. The K F as defined above handles prediction for the S S M with C ^ 0. After all, the S S M can be transformed by modifying the state vector so that C = 0. The case where C is very large can be handled with the K F . However numerical instability becomes a feature. A more efficient and direct approach appropriate to these cases is discussed in the following section. 2.3.5 The smoothing filter Smoothing refers to predicting the state vector a,-, i = 1,..., n, using the entire observa-tion vector y, where y <— SSM. 35 Definition 2.4 Suppose that y <— SSM, and C = 0. Further, assume the KF is run for i — 1 , . . . , n. Then, the smoothing filter is the recursion for i = n,n — where L{ = Ti — K{Z{ and rn = 0, Un = 0. The importance of the smoothing filter derives from the next result and its conse-quences. Theorem 2.8 Suppose that y <— SSM and suppose that i,j = 1 ,2 , . . . m + 1, where m < n . Then, P{cti\y) = Plenty!,... ,ym) + Mi(Lm . . . Li)'rm, Mse(ai,aj\y) = Mse(ai,aj\y1,... ,ym) - cr2Mi(Lm . . . Li)'Um(Lm •••Lj)'Mj, where (Lm . . . Li) = I, if m < i. Fixed interval smoothing The next theorem deals with fixed interval smoothing, where P(cti\y) is required for i = 1 , . . . , n. It gives an efficient algorithm derived from Theorem 2.8. Theorem 2.9 Suppose that y <— SSM, and C = 0. Then, for i = 1 , . . . , n, ri-i = Z'iDi1ei + L\ri, Ui-r = ZlD-'Zi + L'iUiLi, (2.11) P(ai\y) = ai + Pifi-1, Mse\ca\y) = a2(Pi - PU^Pi) (2.12) while for i, j = 1 ,2 , . . . , n + 1, with i < j , Mse(cti, aj\y) = ^P^Lj-x... Li)'(I - Uj^Pj). (2.13) 36 Thus, for fixed interval smoothing, the Kalman filter is used to compute the at-, P, and L{. This "forward" pass is followed by a "backward" or "smoothing" pass to arrive at the P(di\y) and Mse(a,-|y). Fixed point smoothing The fixed point smoothing algorithm computes P(a t|y) for fixed point i and increasing sample size n. A n algorithm is obtained by taking m — n — 1 and i = j in Theorem 2.8: P(ai\y) = P(ai\y1,...,yn.1) + Pi{Ln-1...Li)'Z'nD-1en, Mse(ai\y) = Mse(a , -|yi , . . . , y „ _ i ) - a2Pi(Ln^ ... L^'Z[D~X Z{(Ln-X . ..Li)P{. The implementation requires the concurrent calculation of en, Dn, and Ln, and the accumulation of the factor ( L n _ i . . . Interpolation Interpolation is the procedure where y,- is predicted using y* = ( y i ; . . . ; y , - i ; yi+i; • • •; yn)-Interpolation arises when missing observations are to be predicted, and for cross-validation, for example. Theorem 2.10 Suppose that y <- SSM, and C = 0 and the KF followed by the smooth-ing filter is applied. For i- 1 , . . . , n, define rii = D^ei-K-ri, (2 14) Ni = D~l + KlUiKi. y ' ' 37 Put yu as component I of yi, yl = (yi,;...; j / j _ x ; yi+i;...; yn) and yh as yi excluding component yu. Then, P(yi\y{) = yi-N-'m, Mse(yi\yi) = a2Nr\ P(yu\yli) = yu-N^nu, Mse(yii\yli) = a2N^\ K ' where nu and Nu are component I and Uh diagonal entry of rii and Ni, respectively. Forecasting with the S S M Forecasting can be easily implemented by using the K F and the iterated prediction results. The A;-step ahead predictor of the state, given Yn = (y\,... ,yn), in light of Theorem 2.4, can be iteratively computed using p(an+k\Yn) = p(PK+fc|yn+fc_!)|rn), = Wn+k-id + Tn+k-iP(ctn+k-i\Yn), with Mse given by o-~2Mse(an+k\Yn) = r n + fc_ iMse(a n + A ; _i|y r n)^ + A . _ 1 + i J n + f c - i i J ^ + f c _ i , where P(an\Yn) = an and Mse(an\Yn) - Pn. Therefore, the A;-step ahead observation predictor and its Mse matrix are P(yn+k\Yn) = Xn+k(3 + Zn+kP(an+k\Yn), Mse(yn+k\Yn) = Zn+kMse(an+k\Yn)Z'n+k + a2Gn+kG'n+k. 38 2.4 The diffuse Kalman filter In the previous section, the estimation and prediction tools were given for the case when C in Definition 2.1 is zero. In this section the case where C ^ 0 will be considered. The diffuse Kalman filter ( D K F ) is the tool devised by de Jong (1991) to deal with the S S M when 7 in Definition 2.1 is diffuse. A random vector is said to be diffuse, denoted vectors arise in two ways. First, to model parameter uncertainty and second, diffuseness arises when a nonstationary model is assumed to have applied since time immemorial. The concept of processes which have applied since time immemorial will be defined and treated in the next section. The D K F is a generalisation of the K F and it deals with the estimation of the diffuse parameters of a S S M at the same time that it deals with prediction of the state. The D K F turns two vector recursions in the K F into matrix recursions. D e f i n i t i o n 2.5 The DKF is the KF with the equations for et- and a,-+i, respectively re-C 0 0 , if c*-1 converges to a zero matrix in the Euclidean norm. Diffuse random (2.16) with starting condition A\ = WQ(—B,b). Also, the following recursion is added: Qi+i = Qi + E'i + D~lEi where Q\ = 0. 39 From the previous definition, the D K F is the algorithm: Ei = (XiB^i-X^-ZiAi, Di = Z{P{Z[ + G{G\, Ki = (TiPiZ'i + HiG^D-1, Ai+1 = Wi(-B,b) + TiAi + KiEi, Pi+1 = {Ti-KiZjPiTi + iHi-KiGiW, Qi+i — Qi + E\D~^Ei, i = 1 , . . . , n , with starting conditions A\ = W0(—B, b), Pi = HQH'q, Q\ = 0. The next result shows how the D K F generalises the K F . T h e o r e m 2.11 (de Jong, 1991) Suppose y <— SSM, 7 is fixed (i.e. C = 0), y is normally distributed and the DKF is applied. Let Qn+i = {(S, s); (s', q)}. If S is non-singular then the mle's of 7 and a2 are 7 = S~xx and a2 = (q — s'S~1s)/y#, respectively. The (7, a2)-maximised log-likelihood is y#{ ln («7 2 )} + X>|A| t=i (2.18) 2.4.1 The diffuse likelihood If 7 is random, one has the following result. T h e o r e m 2.12 (de Jong, 1991) Suppose y <— SSM and (7;?/) is normally distributed. Further suppose the DKF is applied and S is nonsingular. Then the mle of c is 7 = S_1x with covariance matrix c r 2 , ? - 1 . The mle of a2 is a2 = (q — s'S~ls)/y*,while the mle of C is zero. The log-likelihood \(y) maximised with respect to c, C, and a2 is - 1 / 2 [y#{ln(o- 2 )}-f-Er = 1 ln|A The next theorem considers the behaviour of A(y) as C — > 0 0 . 40 Theorem 2.13 (de Jong, 1991) Suppose that y <— SSM and ("f,y) is normally dis-tributed. Further, suppose the DKF is applied and S is nonsingular. Then as C —> oo, A(y) + 1/2log \C\ converges to y* log(a2) + log\S\ + Y > g |A| + (q- s'S^s)/^ i=\ (2.19) Moreover, (2.19) is a log-likelihood based on Ny where N has rank y * — 7*, Cov(Ny,*y) — 0 and log\Var(Ny)\ equals the first three terms in (2.19). The log-likelihood (2.19) max-imised with respect to a2 is y*log(a 2) + log|S| + y > g | A | (2.20) 2.4.2 Diffuse prediction The D K F can be used to compute diffuse predictors of a:,- and y,- in the context of the S S M . These are predictors constructed under the assumption that y <— SSM with C —* oo. For i — 1, + 1, let 7i = P ( 7 | y i , . . . , y i _ i ) , at- = P(a,-|yi, • • •, yi-i) and y; P(yi\yi, • • • , y i - i ) - In the next theorem, let Qi = Si Theorem 2.14 (de Jong, 1991) Suppose that y <— SSM where C ^ 0 and the DKF is applied. Let A ; 7 and denote all but the last columns of Ai and Ei. Then, 7 i = (Si + C - ^ i C ^ c + Si), Msetfi) = a^Si + C-1)-1, &i = A t ( - 7 ; ; 1 ) , Mse(ai) = cr2Pi + A^Mse^A'^, yt-Vi = Ei(-v;l), Mse(yi) = a2D{ + E^Mse^E'^. Theorem 2.15 (de Jong, 1991) Suppose that y <— SSM and i — 1 , . . . , n + 1. If Si is nonsingular, then as C —> oo, 7,- —> S^Si, Msefy) —• a^S^1. If the rows of A^ are in the row space of Si, then as C —• 0 0 , on —> A;( —SfS{', 1) and Mse(6ci) —• 41 a2(Pi + A^Si A'iy). If the rows of are in the row space of Si, then as C —> oo, V i {Xt{-B, b) + ZiAi) (-S-Si-1), and Mse(yi) -+ a\Di + E ^ E ' ^ ) . 2.4.3 Diffuse smoothing Smoothing refers to predicting the state vector a;, i — 1 , . . . , n, using the entire observa-tion vector y, where y <— SSM. Smoothing is based on the backwards recursion Mi-r = Z'iD^Ei + LM, (2.21) Ri-! = ZlDr1 Zi + L'iRtLi, (2.22) where A/"n = 0, Rn = 0 and for i = 1 , . . . , n , = T; — -KjZT,, and all other quantities are defined as in the D K F . T h e o r e m 2.16 (de Jong, 1991) Suppose y «— SSM, where C / 0 and the DKF is applied followed by the recursions (2.21) and (2.22). For i — 1 , . . . , n + 1, suppose &i = P(cti\y). then, &i = (At- + PiA^_ 1 )(- 7 t - + 1 ; l) , Mse{ai) = a2(P{- PiRi^Pj+ATi^Msetfn+JjVU^ where j\fi-\fl denotes all but the last column of Ai + PiAfi-i and 7 n + i and M s e ( 7 „ + i ) are as given in Theorem 2.14 with i = n + 1. Furthermore, as C —> oo and provided A / i _ i ) 7 is m £/ie row space of Sn+i, &i and Mse(ati) converge to the previous expressions with %+i replaced by S ^ + i S n + i and Mse(%+i) replaced by cr2S~+1. 42 2.5 State space models which have applied since time immemorial In this section it will be assumed that the S S M (2.1) is index invariant, that is, for i — 1 , . . . , n, Wi = W, Ti — T and Hi = H, where W, T, and H do not depend on i. Definition 2.6 A index invariant SSM is said to have applied since time immemorial if the state transition equation = W/3 + Tcti + Hui is assumed to hold for i = 0 , - 1 , - 2 , . . . . Thus, under the time immemorial assumption, the initial suppositions about ctrj are superseded. In what follows, a root of a matrix T is an eigenvalue of T. A root, A, of T is called stationary if |A| < 1, and non-stationary otherwise. Theorem 2.17 (de Jong and Chu-Chun-Lin, 1994) Suppose a SSM is assumed to have applied since time immemorial. Further, suppose U is any matrix such that UTU'1 = diag(P, Q) where P has all non-stationary roots and Q all stationary roots. Suppose S is the last set of rows in U with the number of rows equal to that in Q and put U = (U\, L ^ ) - 1 where U2 has the same number of columns as Q. Then for t = 1,0, — 1,.. [Var(ai)r = a^S'M'S = a'2 l im UUxUi + U2MU')~ , K—•OO where M satisfies M = QMQ' + SHH'S'. Consequences of the previous theorem worth noting are the following. 43 • If T has all stationary roots, then U\ is null, C/2 = S = I and Q = T. Hence, if T has all stationary roots, Var(a{) = a2M, where M satisfies the matrix equation M = TMV + EE'. This result is also given in Anderson and Moore (1979). • If all the roots of T are non-stationary then U2 is null, U\ — I and [Var(a: t-)]_ = 0. Stating things loosely, one could say that 7ar(a , ) = <T2KI, where K — > 0 0 . This result is often arrived at intuitively by noting that in the non-stationary case Var(a,-) is arbitrarily large. • If T has both stationary and non-stationary roots, as is usually the case in empirical non-stationary cases, the time immemorial argument leads to an arbitrarily large Var(a,-). However the divergence is specialized in that [Var(a;)]~ is a finite matrix not necessarily equal to zero. • A matrix U as specified always exists. This is a consequence of the Jordan decom-position discussed for example in Shilov (1971). 2.5.1 Prediction and likelihood evaluation for non-stationary models The problem of initialising the D K F , that is the problem of modeling the initial state ao is addressed in this section. This is material from de Jong and Chu-Chun-Lin (1994). D e f i n i t i o n 2.7 The diffuse SSM is the SSM as in Definition 2.1 except that ao — a + AS + f / _ i u _ i and (3 = b + B6 where (a; b) is fixed and known. (A; B) is of full rank, E(6) = c, Var(6) = a2C with C non-singular unless C — 0 and 8 is uncorrelated with ( u _ i ; u 0 ; u i ; . . . ; u „ ) . 44 The vector 8 models unknown aspects of the starting state ao and the parameter vector 0. The cases of interest are where C = 0 and C~l = 0. The case C = 0 models a fixed but unknown 8. The case C _ 1 = 0, denoted C —• oo, indicates that (cvo; 3) is random with an arbitrary large, albeit specialized, covariance matrix. In other words (cto; 3) is said to be diffuse. If (A; B) has full row rank then (ao-, 3) is fully diffuse; it is partly diffuse otherwise. T h e o r e m 2.18 (de Jong and Chu-Chun-Lin, 1994) Suppose that yi is generated by a time immemorial SSM with 8 known and assume the notation of Theorem 2.17. If RW0 = 0 then yi is generated by the diffuse SSM with a — U2{I — Q)'1 SW 3, A = U\, B = 0, #_i - U2Mil\ In the case of unknown 3 the diffuse aspects of the model can be structured in the following manner. 1. Using the time immemorial argument, the initial state cvo is assumed to have the specification a0 = U2{I - QYXSWQ + + U2M1/2u-U where 8\ ~ (0,a2kl) with k —• oo. 2. The parameter vector 0 is structured as 3 = b + B282 where 82 ~ (0, a2kl) is the diffuse vector modeling unknown 3. The vector b and matrix B2 capture linear restrictions on 3. 3. The specifications in 1. and 2. are combined to yield, a0 = U2(I-Q)-1SW(b + B282)-rU181-rU2M^2u.1 = U2{I - Q^SWb + {Ux, U2(I - Q)-1SWB2}6 + U2M1/2u.u 45 and hence (A,a) = {U1,U2{I-Q)-1SW(B-2,b)}, H.1 = U2M1^. (2.23) If W f3 is assumed to be zero in 1., the time immemorial argument is used to generate an appropriate Var(ao). In this case, (A,a) = (U-\,,0). The next result shows how the D K F is initialised in the non-stationary case and solves the problem of the D K F initialisation. Theorem 2.19 (de Jong and Chu-Chun-Lin, 1994) Suppose that is generated by a time immemorial SSM. Suppose (Ui,U2) and M are derived from T and H as described in Theorem 2.17. Then the DKF is initialised with Ax = W(0, B2, b) + T(A, a), Pt = TU2MU'2T' + RE\ Qx = 0 where (A, a) is as specified in (2.23) or (A, a) = (/7i,0). The DKF equations for E{ and Ai+1 of the DKF (2.17) become Et = (0,-XiB^yi-Xib)- ZiAi Ai+1 = (OWBitWity + TiAi + KiEi. Note that if B2 is null and T has all stationary roots then A\ = Wb + Ta is a vector and P i = M. In other words, the D K F reduces to the ordinary Kalman filter. 2.6 Cross-validation and related methods of estima-tion A cross-validation criterion can be used for statistical inference. In its original form, it consists of the division of the data sample (of size n) into a "construction" subsample (of 46 size n — 1) and a "validation" subsample (of size 1) in all, re, possible ways. The idea is to optimize on the "construction" sample by testing on the "validation" sample, repeating this for every possible case. Because of its "distribution free" nature, cross-validation is preferred as a method of estimation to maximum likelihood when the Gaussianity of the errors in the S S M is believed not to hold. Suppose one has re pairs of data values S = {(xi,y(xi)), i = l , . . . , re } . Let S-i = S \ {(xi,y(xi))}, i = 1 , . . . ,re, and y(x,9, S-i) be the predictor of y given x, 9 and S-i. Here we are assuming that the predictor of y depends not only on x and the observed sample but also on a (generally unknown) parameter 9, where 9 £ A C Rk, k = 1,2,.... Define the cross-validation criterion as n CV{9) = ^ 2(y(xi) - y(xi, 8, <?_,))2. (2.24) i = l Then, an estimate of 9 will be 8 such that 8 minimizes CV{8). The estimate 8 is the value of 8 under which the missing data points yi are best predicted by the remainder of the data. 2.6.1 Cross-validation for state space models In the context of state space models, predictors will be best linear predictors as defined in Definition 2.3. 47 Theorem 2.20 Suppose that y <— SSM, and C = 0, and the KF followed by the smooth-ing filter is applied. Then the cross-validation criterion (2.24) 45 - PfaWiVi - PW) = It^ryNr^, (2.25) where nt- and Ni are as in Theorem 2.10. Proof . The result follows immediately from Theorem 2.10. • 2.6.2 Generalised cross-validation and state space models Suppose that the signal plus noise model holds, Y = S + 6, where Var(Y) = E y , Var(e) = a2I, Cov(s,e) = 0, and Cov{s,Y) = E s y . Then P(s\Y) = HY, where H, the so called hat matrix, is H = E ^ E " 1 . The generalised cross-validation criterion is defined as GCV = -{Y - P(s\Y))\Y - P(s\Y))/ [1 - t r (#) /n] 2 . (2.26) n Note that one can rewrite (2.26) as GCV = -Y'{I - H)\I - H)Y/ [1 - t r (#) /n] 2 . (2.27) n From the fact that, according to the specification of the signal plus noise model, T,sy = Var(s) and E y = Var(s) + cr 2/, one has = ( E y — E S 2 / ) E 2 / 1 V 48 Also, tr(J -H)=n- tv(H) = n ( l - tv(H)/n). Therefore, the generalised cross-validation criterion can be written as GCV = nY' [Var{Y)}-2 Y/ [tr ( ( y a r ( F ) ) " 1 ) ] 2 . (2.28) From de Jong (1988) one can also write the generalised cross validation criterion as EILi D [(yu - P(yii\yli))/Var(yii - P(yii\yu))]2 E ? = i E / l / V a r ( y K - P ( y K | y « ) ) ] 2 and so G C V can be computed with a pass of the Kalman filter followed by the smoothing filter as GCV = X > ; ) ' n ; t'=l 1 2 £ t r ( J V t ) (2.30) where the quantities n4- and Ni were defined in Theorem 2.10. 2.6.3 Generalised cross-validation and diffuse state space mod-els If there are diffuse components in the S S M , then the generalised cross-validation criterion can also be computed using the D K F and the diffuse smoothing filter, as it is shown in de Jong (1988). Let us assume now that the diffuse S S M , as per Definition 2.7, holds. The following equation should be added in the diffuse smoothing filter. Mi=Dr1Ei-K'iAfi, where jV% is given in (2.21). It is shown in de Jong (1988) that as C —> oo, the G C V criterion converges to GCV = £ ( M t ( - 7 n + 1 ; 1))' (Mt(-%+l] 1)) j jx>(JVi + M i 7 5 " + 1 M ^ ) | (2.31) where TV,- is given in Theorem 2.10 and M ; 7 denotes all but the last column of M,-. 49 2.7 Summary This chapter consisted of the review of existing material on discrete state space models. The reason for this review is that this material will be crucial for the understanding and development of next chapter. In this chapter the state space model has been defined. The intuitive ideas behind such stochastic constructs have been illustrated. It has been illustrated how starting from a global functional form for the signal one can, by the appropriate placement of disturbances, construct local models that will reflect local characteristics of the data. The estimation and prediction principles and tools to be used with the state space model were presented. Also, a generalisation of the state space model, the diffuse state space model was reviewed. The corresponding generalisations of the prediction and estimation tools were also given. Finally an alternative to the maximum likelihood estimation procedure, the cross-validation and generalised cross-validation techniques were introduced. 50 Chapter 3 Continuous state space models Although most variables are measured at discrete intervals, there are many advantages of formulating a continuous model for smoothing. A model formulated in continuous form does not depend on any intervals at which observations are made. Once data are collected, from the continuous model one may derive a discrete model to suit the actual observations of the variables. Some economics researchers, such as Bergstrom (1976, 1984), argue that, in fact, some real processes generating observations are continuous hence turning continuous models into more fundamental than discrete ones. Usually they refer to continuity of processes in time. In this thesis, we do not want to argue about the physical fulfillment of continuous models but on the statistical arguments and benefits of their use. It needs to be stressed again that the fact that a model is stated in continuous form does not imply the process in question happens in continuous form. It is the view of the author that a mathemat-ical (stochastic) model has the purpose not to exactly dictate how reality occurs but to understand its manifestations. In this sense, a stochastic model will be used as a mere tool to extract as much information as possible from collected data in order to better 51 comprehend the processes that originated them. In this chapter continuous state space models will be introduced as a tool for smooth-ing. A continuous state space model is a continuous signal plus noise model where the signal is a continuous stochastic process that satisfies a stochastic differential equation. Section 3.1 contains a formal definition of continuous S S M . Continuous state space models are conceived using the same basic idea that discrete state space models are built with. It is assumed first that the signal satisfies a general functional form, in this case given through a differential equation. Then, by an appropriate placement of disturbances, one constructs a more flexible stochastic model for the signal. Sections 3.2 and 3.3 indicate how to implement estimation and prediction procedures in the context of a continuous S S M . These two sections are a new contribution, for the results therein can be used for any continuous state space model. So far, only the case where the process defined by the state equation is stationary has been dealt with in the literature. Section 3.4 deals with L-splines, a smoothing methodology that, although having origins in the area of nonparametric statistics, can be cast into a continuous state space model. This is not surprising since the idea behind the L-spline technique is that the signal has a functional form that satisfies a differential equation. Departures from this functional form are allowed since the signal is estimated by minimising a penalised least squares criterion. The penalty of this least squares criterion is the integral of the differential operator that the idealised initial form of the signal satisfies. So, big departures from this initial form of the signal are penalised. A n example is presented that compares, using real data, the approaches of L-splines and continuous state space models in data analysis. 52 Section 3.5 reviews and expands the topic of continuous structural models. A n example is analysed to show how an L-spline can be cast into a continuous structural model. A comparison with another nonparametric methodology (polynomial smoothing spline) is offered and the advantages of using a stochastic model for data analysis and comparison of models are illustrated. 3.1 Definition of continuous state space models We will consider continuous state space models where the state equation is a stochastic differential equation for the state vector, namely • e(x) is a white noise process with variance a2. • Wx is a standard Wiener process. For a definition of white noise and Wiener process see Chapter 4. Some basic theory of stochastic differential equations will be essential to treat continuous state space models. For some background and references on stochastic differential equations, see Appendix A . 3.2 Discrete version of a continuous state space model Suppose that one has the model (3.1). Then, from Theorem A.3 , for 8 > 0, X(x)0 + Z{x)a(x) + G{x) e(x) Aa(x) dx + H dWx x > Xi, (3.1) where (3.2) 53 Hence, for the data points (a;,-, y t ) , i — 1 , . . . , n, one has the following discrete state space model f yi = Xi(3 + ZiOti + Giti . > \ a i + 1 = TiOLi + Ui where yt - y(xi), Xi = X(xi), Gi = G(xi), ai = a(xi), Zi — Z(xi), i = 1 , . . . , re, T- = ehiA, with hi — x,+1 — Xi, i = 1,..., re — 1. The u,'s are independent random variables with m = fh' e(h>-t)AHdWt. (3.4) Jo From (3.4) and Theorem A.1 it follows that E(Ui) = 0, and Var(Ui) - e{h'-s)AHH'^hi-^A' ds. (3.5) Jo Comparing model (3.3) with (2.1) one notes that the term Wi/3 in the state equation is not present. This is due to the formulation of the stochastic differential equation defining the state. Hence, we assume here Wi in (2.1) to be zero. 3.3 Estimation using the continuous state space model In order to implement the D K F , the tool that will be used for estimation, one needs to compute the state transition matrix, the state error variance and deal with initialisation issues. 3.3.1 The state transition matrix To implement the D K F , one needs to compute the state transition matrix, Tj- = eh'A, given the matrix A. 54 The definition, properties and examples of matrix exponentials are given in Appendix D . In most of the examples presented in this thesis, the matrix exponential will be computed explicitly. There are ways of computing the exponential of a matrix numerically. See Appendix D for references. 3.3.2 Variance matrix of the state error In order to implement the algorithms, one has to compute Var(ui), as in (3.5). This is a more difficult task since it involves integration of the matrix exponential. It is known (see e.g. Arnold, 1974 or Theorem A.5 in Appendix A) that if the eigen-values of the matrix A in (3.1) all have strictly negative real parts then the stochastic process defined by the state equation da(x) = Aa(x) dx + H dWx is stationary with mean zero and variance K, where K is the solution of the matrix equation AK + KA' = —BB'. This result is appealing for it avoids integrating the matrix exponential to compute the state error variance. Instead a linear matrix system has to be solved. The case where some or all the eigenvalues of A have real part greater than or equal to zero has not yet been dealt with, at least in the reviewed literature. This case is of special importance since for most continuous structural models and continuous S S M for L-splines (see Sections 3.4 and 3.5) the matrix A has some or all of its eigenvalues with this property. In this section we develop a new result on how to compute Var(ui) by solving a linear matrix system. It is shown that all that is needed is to first compute the matrix exponential of A and then solve a linear matrix system. L e m m a 1 Let (3.6) 55 Then, V(h) satisfies V{h)A' + AV{h) = ehAHH'ehA' - EE'. (3.7) P r o o f . Note that by a change of variables in (3.6) one obtains rh V(h) = / JO esAEE'esA> ds. Also, one can write V(h) = ehA \^j\-sAEE'e-sA' Hence, from (3.8) and (3.9), £ esAEE'esA> ds = ehA Ijf ds esAHH,e-sA< d s ,hA' (3.8) (3.9) (3.10) Differentiating with respect to h on both sides of (3.10) one obtains the desired result. • R e m a r k 3.1 It is known that the matrix equation BX + XC = D (3.11) has a unique solution, namely /•oo X = - etBDetCdt, (3.12) Jo provided that the integral on the right hand side of (3.12) exists (see e.g. Bellman, 1960). A few noteworthy consequences are the following. 56 If (3.11) has a solution X and D is symmetric, then from (3.12) it follows that if C = B then X must be symmetric. If either B or C have all stationary roots, then the solution to (3.11) exists and is unique. In order to prove this claim, using Theorem D . l in the Appendix, it will be enough to consider that B and C are elementary Jordan blocks (see ( C l ) , in the Appendix). Hence, etB = eXtM(t) and etC = e" 4 M(t) , where M(t) is a matrix whose entries are either constants or multiples of tk, for some positive integer k (see Appendix D) . Suppose that Re(A) < 0 and Re(/j) > 0. Then, Re(A - fi) < 0 and I [°° etBDetC dt\ < re^-^ \M(t)DM(t)\ dt < oo, \Jo I Jo where \A\ denotes the matrix of absolute values of the entries of the matrix A. If B = A and C = A', for some matrix A, then one can transform (3.11) to an equivalent equation. Let P be the matrix such that PAP'1 = Si ': 0 (3.13) o ; s2 where Si has all non-stationary roots and S2 has all stationary roots. Such matrices P, Si and 62 exist and are given by the Jordan's canonical form of the matrix A (see e.g. Shilov, 1971). Hence, if B = A and C = A'', (3.11) is equivalent to A*X* + X*(A*)' = C*, (3.14) with A* = PAP'1, X* = PXP' and C* = PGP'. 57 Partition X* and C* conformaily to diag(5'i, 2^), such that X* x*u X12 X21 • X22 c* = C\2 C21 ^22 Then, (3.14) is equivalent to the system of equations S\X\X + x\x s[ S\Xl2 + x^s'2 S2X21 + X21 S[ S2X22 "F X22S2 U12) 2^1) 2^2-(3.15) (3.16) (3.17) (3.18) Note that (3.16) and (3.17) are equivalent and that (3.16), (3.17) and (3.18) always have a unique solution since Si has all stationary roots. • In reference to Lemma 1, write V(h) = Vu(h) V12(h) v12(hy ; v22(h) where the partition of V(h) is conformal to (3.13). Then, V\2(h) and V22{h) can always be uniquely determined. V\\(h) may or may not be determined. This , as we shall see, is immaterial in terms of the D K F , under the time immemorial assumption. First note that (3.1) is equivalent to y(x) = X(x)0 + Z*(x)ct*(x) + G(x)e(x) da*(x) = (PAP-1)a*(x) + H*dWx where P is as in (3.13), Z*(x) = Z{x)P~1: a*(x) = Pa(x) and H* = PH. 58 Let x be fixed and h > 0. Then, a*{x + h) = d^g(ehS\ehS2)a*(x) + u(h), and so Var(a(x + h)) where ehSlNiiehS[ : ehS2N2lehSi : Var(a(x)) = ehSl Ni2ehiii hs< ehS2N22ehS2 + Vuih) :• v12(h) V12(h)' \ V22(h) N-11 12 N21 : N22 Now, ehSlNuehS'i and Vu(h) diverge to infinity in Euclidean norm as h —>• oo since S-L has all stationary roots, whereas ehS2N22e,hS'? —> 0, ehSl N12ehS2 —• 0 as h - * oo. This can be easily seen by first assuming that S\ and 52 are elementary Jordan blocks and then extending the result using Jordan's canonical form of a matrix. Hence, under the time immemorial assumption, it is immaterial how one determines Vn(h). a Remark 3.2 It is a consequence of the proof of the previous theorem that V(h) satisfies the determin-istic differential equation dV(h) dh = V(h)A' + AV{h) + EE'. (3.19) 59 Applying the vec operator (see Appendix) on both sides of (3.19) one obtains dvec(V(h)) dh = Bvec(V(h)) + vec(HH'), (3.20) which is a deterministic first order linear equation for vec(V(h)) with initial condition V(0) = 0 and B is defined as B — (A <g> 7) + (7 ® A). The solution to this differential equation is vec(V(h))= [k etB dtvec(HH'), (3.21) Jo which in the case when A is invertible (and hence B is invertible) vec(V(h)) = B'1 (ehB - i) vec(7777'). Now, observe that from the properties of the Kronecker product (see Appendix) ehA®ehA = (I + hA + h2/2\A2 + ---)®{I + hA + h2/2\A2 + ---) = I + hB + h2/2\B2 + ---= ehB. • In the case where B is not invertible one can use the following result. T h e o r e m 3.1 The matrix V(h) defined in (3.6) satisfies Bvec(V(h)) = [(T(h) <g> T(h) - I)]vec(HH'), (3.22) where B = (A® I) + (I ® A). 60 P r o o f . From Lemma 1 one has that V(h) satisfies V(h)A' + AV(h) + HH' = ehAHH'ehA' (3.23) Applying the vec operator (see Appendix) on both sides of (3.23) one obtains Bvec{V(h)) = (ehA <g> ehA - I) vec{HH'), (3.24) which is what we wanted to prove. • The importance of the previous theorem is that it implies that all that one needs to compute V(h) is first to compute T(h), i.e. the exponential of the matrix hA, and then solve the linear system (3.22). In the case when B is not invertible, there will be more than one solution to (3.22) and a positive definite one should be sought. We know already, from Remark 3.1 the nature and source of indeterminacy in V(h). As it was mentioned in this Remark, under the time immemorial assumption, the undetermined components can be arbitrarily chosen so that the matrix V(h) is positive definite. To implement the above computational procedure, one could use some of the many available mathematical software packages. The ideal package should have a function that computes Kronecker products (to obtain the matrix B), matrix exponentials (to obtain T(h)) and solves linear equations. Maple computes matrix exponentials and solves linear equations. It does not compute Kronecker products. Another package worthwhile investigating is Octave. Octave is a high-level language, primarily intended for numerical computations. Its linear algebra section has the three required functions (exponential of a matrix, Kronecker products and solve linear systems). 61 E x a m p l e 3.1 (The variance of the state error for an exponential decay model) For details about the exponential decay model, please refer to Example 3.4. In Exam-ple 3.4 the variance of the state error is computed explicitly. To illustrate Theorem 3.1, we will compute the variance matrix using the method proposed in that theorem. In this case, A Hence, 0 0 0 0 T(h) 0 1 0 ehe B = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 and T(h)®T{h)-I = + 0 0 0 e h e - l 0 0 0 0 and H = T/0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 - 1 1 e w - l 0 0 0 Then, assuming that V(h) = a(h) b(h) b(h) c(h) 0 / o \ b(h) b(h) V 2c(h) J , (3.22) is in this particular case 0 l - e w l - e h e Ve 2h0 or b(h) = £ ( l - e " ) , <h) = g r ( e ^ - l ) . As we can see there is one indeterminacy and a(h) has to be chosen such that V(h) is positive definite. The choice a(h) = (r/6)2h makes V(h) positive definite. Note that V(h) coincides with the expression given in Example 3.4. • 62 idealized past times -< .. *• actual observations 3 1 n i -I i i i i i Xi — 2 X\ — 1 Xi X2 X3 X4X5 Figure 3.1: The time immemorial assumption. 3.3.3 Initialization of the Diffuse Kalman filter when using a continuous state space model In order to carry out estimation within the context of the continuous state space model (3.1) one needs to specify the properties of a(xi) = ax, E(cci) and Var(ai). For this purpose, similar to de Jong and Chu-Chun-Lin (1992), it will be assumed that the process a(x) satisfies the time immemorial assumption, that is, for x — xi — 1, xi — 2, xi — 3 , . . . a(x + 1) = eAa(x) + u(x), where E(u{x)) = 0 and Var(u(x)) = V(l) = ft evAHH'evA' dv. Then, E(a(x + k)) and Var(a(x + k)) will be considered as k —y 0 0 , over the integers, given that the process a(x) starts at x — Xi. This is equivalent to studying E(a(xi)) and Var(a(xi)) under the time immemorial assumption. A l l the results stated in de Jong and Chu-Chun-Lin (1992), summarised in Chapter 2, apply directly to this situation. E x a m p l e 3.2 (mth order linear stochastic differential equations.) Suppose that Y(x) = s(x) + e(x) where e(x) is a white noise process with variance <r2 and 63 the signal, s(x), satisfies the stochastic differential equation $<m>(x) + a^s^-^ix) + •••+ a0s(x) = r u(x) where u(x) is a standard white noise process. Then, from Remark A . l , one can write the following continuous state space model for Y Y{x) = ( l , 0 , . . . , 0 ) a ( s ) + c(x) da(x) = Aa(x)dx + (0,...,0,T)dWx where ct(x) = (s(x), a ( 1>(x),..., s( m _ 1 )(x)) and A 0 1 0 0 0 0 1 0 0 0 —a0 —ax 1 0 0 1 (3.25) (3.26) Using the same idea with which structural models in the discrete case were built, i.e. placing disturbances appropriately in order to obtain a more flexible model, one could think of the following model, instead of (3.25) Y(x) = ( l , 0 , . . . , 0 ) a ( x ) + e(x), da(x) — Aa(x)dx + H dWx, (3.27) where H = d i a g ( r i , . . . , r m ) , T,- > 0, i = 1 , . . . , m. In this way we can put a mth order stochastic linear equation as a first order system. • Examp le 3.3 (A fourth order linear stochastic differential equation to model a signal with the "linear trend + cycle" structure.) 64 In this example the initialisation issues will be dealt with in a specific case. Suppose that, in the previous example, s(x) satisfies the stochastic differential equation S(iv\x)+u2sM(x) = u(x), that is, m = A, ao = ai — a3 = 0 and a2 = —u2. To simplify computations it is convenient to transform the state equation thereby clearly revealing the root structure of the state transition matrix. If the matrix P is defined as P = one has that P - 1 = and P-XAP = 1 0 --1/u3 0 0 1 0 -1/u2 0 0 1/u 0 0 0 0 1 r 1 0 1/u2 0 1 0 1 0 1/u2 0 0 U 0 . o 0 0 1 ' 0 1 0 0 " 0 0 0 0 0 0 0 u . 0 0 -u 0 . Hence, model (3.25) is equivalent to r Y(X) = \ ala*(x) = ( l , 0 , - ^ , 0 ) a * ( x ) + e(x) A*a(x) dx + H* dWx where a*(x) = P~ 1 a(x ) and H* = P-XH, that is H* T l 0 ^ 0 0 T 2 0 % 0 0 UT3 0 0 0 0 r 4 j (3.28) 65 One also has that " 1 t 0 0 " 0 1 0 0 0 0 cos(tu>) sin(iu>) . 0 0 — sin(£o;) cos(tu;) H*(H*Y 0 UJ 0 0 rZ + & 0 0 0 0 0 etA*H*(H*yet(A*)> = Mt(t) M3(t)' where M3(t) M2(t) Mx(t) ( r | + ^ + ( r a 2 + ( r 2 2 + Hence, with [SetA'H*(H*)'et(A'ydt = Jo Vi(6) = Vi(S) V3{8) (V3(8))' V2(8) V2{8) ^^Sin(2u8) + l(^ri + r4)8 sin2(a;<5) 2w (T?-W2T?) . 2 / c\ 4 4 c s in(2a;*) + ^ 66 and + sinM) - %6 cosM) + cosM) - sinM) - © + 2 2 -^(cosM) - 1) & sin(w5). Note that since A * = d i a g ( M i , M 2 ) , where 0 1 0 0 and M 2 = 0 u -u 0 eA* has all non-stationary roots. Therefore the matrix U of Theorem 2.17 is U the D K F should be initialised, according to Theorem 2.19, with Ax = (e^ ,0) , P a = H*(H*)', and Qx = 0. I and • 3.4 Smoothing L-Splines There is a continuous S S M , such that the solution to the signal extraction problem is the smoothing L-spline, as will be defined below. This continuous S S M has the feature that it always assumes completely diffuse initial conditions for the state vector, regardless of the state transition matrix. As has been discussed, the initial conditions are not always diffuse and they can depend on the state transition matrix. Fixing initial conditions to be completely diffuse regardless of state transition matrix will be shown to give inadequate results in some cases. In this section, L-splines will be reviewed. The attention will be put on published work that attempts to link the functional analytic side of L-splines with stochastic processes, particularly with state space models. For a good source of reference on L-splines see, for example, Weinert, Byrd and Sidhu (1980), Wahba (1990), Ansley, Kohn and Wong (1993) and the references therein. Recent 67 work by Ramsay (1996b) also addresses L-splines, although from a non-model based ap-proach. Kimeldorf and Wahba (1971), obtained a correspondence between spline smooth-ing and stochastic estimation. This correspondence, however, does not lead to useful computational algorithms. 3.4.1 Definition of L-splines Suppose that a regression function / is observed with noise, so that, if y(x) = f(x) + e(x), where e is a white noise process with variance cr2, only yi — y(xi), i = l , . . . , n , are observed. The interest lies in obtaining an estimate of / . Suppose further that / is similar in shape to a parametric function g(-; 9) which satisfies that Lg — 0 where L is the the rath order differential operator where D(f) = df /dx. It is assumed that the functions a,-(-; 9), i = 1 , . . . , ra, are ra times differentiable. Such parametric functions appear frequently in model growth and decay curves in ecology, product diffusion in marketing, finance problems, etc. To allow departures from the parametric model, it is required that the L-spline solves the following optimization problem for given 9 and A. The space W ^ f O i 1] 1S the rath order Sobolev space and is defined as L = Dm + ai(-; 9)Dm~1 + ••• + a m(- ; 9)1, (3.29) (3.30) W™ = If : is absolutely continuous, j = 0 , . . . , ra — 1, and / W 6 I 2 [ 0 , 1 ] } . 68 If ai = 0, 1 = 1 , . . . , ra, then the solution to (3.30) is a piecewise polynomial of order 2ra — 1 called the polynomial smoothing spline. See Chapter 4. 3.4.2 A stochastic model for L-spline smoothing A S S M for L-spline smoothing was obtained by Weinert, Byrd and Sidhu (1980) which was further studied by Ansley and Kohn (1986). A better formulated model was given by Kohn and Ansley (1987) and further explored by Ansley, K o h n and Wong (1993). This model is based on the following result. T h e o r e m 3.2 (Kohn and Ansley, 1987) Let yi = /(#;) + e,-, i = l,...,n, where f € W 2 m [0 , l ] satisfies Lf(x) - au(x), with L being the differential operator defined in (3.29) and e and u are independent stan-dard white noise process. Also, assume that (/(0), / ' ( 0 ) , . . . , / ( r n _ 1 ) ( 0 ) ) ' is normally dis-tributed with mean 0 and variance hi, k > 0, and that this vector is independent of u and the ti's. Then, the minimizer of (3.30) is given by f(x) = l im E[f{x)\yu...,yn). k—>oo As a consequence of Theorem 3.2, the solution to the signal extraction problem using the following model gives the smoothing L-spline: for x > x\, y(x) = Za(x) + Xcre(x), (3.31) da(x) = A(6,x)a(x) + aHdWx, (3.32) 0 Im-l where Z = c ( l ,0 , . . . ,0 ) , = (0, . . . ,0,1), and A{6,x) = " , with am—l\X) ... uoyx J Im-\ denoting the (m — 1) x (m — 1) identity matrix and 0 the (ra — 1) dimensional vector 69 with all entries null. The initial conditions are a(xi) ~ N(0,klm) with k large. In order to obtain the smoothing L-spline, the initial conditions always have to be distributed N(0,klm), with k large, regardless of the state transition matrix A. This is what is known as completely diffuse initial conditions. If one assumes the stochastic behaviour (3.32) for the signal a then the initial conditions will not always be completely diffuse. Depending on the matrix A, the stochastic process defined by (3.32) could be stationary, non-stationary (giving completely diffuse initial conditions) or partially non-stationary (which gives partially diffuse initial conditions). If one identifies, based on physical reasons or exploration of the data, that (3.32) is adequate to describe the behaviour of the signal, one should stress this choice of a model rather than the desire to solve an optimisation problem. It should be the data dictating the degree of smoothness of a curve rather than the methodology chosen. The following example attempts to show the effects of assuming completely diffuse initial conditions when in fact partially diffuse initial conditions are dictated by the tran-sition equation (3.32). Also, it shows the importance of having a stochastic model for smoothing since this allows the comparison of models. Examp le 3.4 (Modeling exponential decay.) The ideas exposed above are illustrated by considering the example of Ansley, Kohn and Wong (1993, Section 5). The data are from Draper and Smith (1981, p.476) and can also be found in Appendix E of this thesis. The investigation involved a product which must have a fraction 0.50 of available chlorine at the time of manufacture. It is known that the fraction of available chlorine in the product decreases with time. In the eight 70 weeks before the product reaches the consumer a decline to a level 0.49 occurs but many uncontrolled factors then arise (such as warehousing environments, handling facilities). Cartons of the product were analysed over a period to provide the data. There are a total of 44 observations covering weeks 8,10,12, . . . ,42. The available chlorine is modeled as a function of time by an exponentially decaying curve. If we let Y(x) denote available chlorine at time x, then it is proposed to model Y as Y(x) = s(x) + e(x) where e(x) is a white noise process with variance a2 and the signal, s(x), satisfies the stochastic differential equation, s(2)(x) - 9s(x) = TU(X), where 9 < 0, r > 0 and u(x) is a standard white noise process. Then, from Example 3.2, one can write the following continuous state space model for Y, Y(x) = (1 0)a(x) + G{x)e(x) da{x) - Aa(x) dx + (0, r ) ' dWx where a(x) = (s(x),s'(x))' and A = Defining P as P 0 1 0 9 1 1 0 9 one has that p-i 1 -1/9 0 1/9 and P-XAP 0 0 0 9 Hence, model (3.33) is equivalent to Y{x) = (1 l)a*(x) + G(x)t(x) da*(x) = A*a*(x)dx-rT/9(-l,l)'dWx (3.33.) (3.34) where a* = P 1a. 71 Note that 1 0 0 ete and etA'HH\etA*)' = (r/oy -(T/e)2ete -(r/eyete (r/eye™ where H = ( r / 0 ) ( - l , 1)'. Given that the weeks at which data on available chlorine was recorded are equally spaced, one has the following discrete state space representation of (3.34), Yi = (1 l)a* + met - ~A\ (3.35) i t - + i - oti + Ui where Yi, ti, ai denote respectively Y (the average), e and a at week i , n ; is the number of observations made at week i, i — 8 ,10,12, . . . , 40,42 and u^'s are mean zero i . i .d . , independent of e,-'s with Var(ui) (r/ey -T2/e3(ee -1) -T2/93(ee - 1) T2/(263)(e26 - 1) Note that eA* has one non-stationary root and one stationary root. If one were to apply Theorem 2.19 to initialize the D K F then it should be initialized with 1 0 0 0 , and P i (r/ey -r2/e3{e9 -1) -T2/63(ee - 1) -T2/{203) If instead one uses the L-spline approach where the initial state, regardless of the assumed model, is assumed to be fully diffuse, the D K F is initialized with Ai - e A* 1 0 0 0 1 0 and P i = ( T / 0 ) 2 -r2/63{ee - 1) -T2/63(ee - 1) T2/(263){e2e - 1) The cases considered for estimation are: I) The initial conditions are assumed to be fully diffuse (L-spline). 72 II) The initial conditions are assumed to be partially diffuse. In all the figures that follow, confidence intervals are constructed pointwise and con-nected by straight lines. Also, a dot represents a datum, solid lines are predicted signals or innovation residuals and dashed lines are predicted error bands. Following are the results of estimation. Table 3.1 shows the estimated values of the parameters, obtained by maximising the likelihood. As one can see, the estimates of the parameters are rather different in the two considered cases. The assumption of fully diffuse initial conditions sets the value of the state error to zero and of the smoothing parameter (the ratio of the measurement error variance to the state error variance) to infinity. This is not the case when partially diffuse initial conditions are assumed in which case a finite smoothing parameter and positive measurement and state errors are obtained. Case Estimated parameters a2 T A = a2IT2 6 I 2.8 x l f r 5 0 oo -0.161 II 3.7 x I O " 5 2.4 x I O " 6 15.32 -0.024 Table 3.1: Estimated parameters of chlorine data model for two different modeling situ-ations. Figure 3.2 shows the predicted signals together with the pointwise confidence intervals joined by straight lines. Figure 3.4 shows the two predicted signals superimposed. B y examining these pictures one sees that the predicted signals are not very different to the naked eye and that the confidence bands for Case I (L-splines) are slightly narrower than for Case II. These confidence bands do not seem to be dictated by the data since they leave too many points outside of the area enclosed by them. In Figure 3.3 plots of the series of innovation residuals are shown, for each of the considered cases. From this figures 73 Case I Case II Figure 3.2: Predicted signals with 95% confidence bands of chlorine data using two dif-ferent models. one can see that the residual at week 18 is unusually large. Also, there seems to be more uncertainty in the predictions of the innovations in Case I than in Case II. The conclusion is that assuming fully diffuse initial conditions, when in fact the choice of the model dictates that they are not, has some effect in the estimation. If partially diffuse initial conditions asre used, then more data driven confidence bands for the signal and less uncertainty in the innovations are obtained. Note how the use of diagnostic tools permitted a simple and clear comparison of the models considered in this example. 74 Case I Case II 15 20 25 30 35 40 45 15 20 25 30 35 40 45 Week Week Figure 3.3: Predicted innovation residuals with 95% confidence bands for chlorine data using two different models. 10 20 30 40 Week Figure 3.4: Predicted signals for Case I and II, superimposed. 75 3.5 Some continuous structural models Motivated in economic applications, Harvey and Stock (1986) discuss structural time series models by continuous time. A signal plus noise model of the form y(x) = s(x) + ae(x) is assumed, where e(x) is a standard white noise process. For the linear trend component, if one assumes that the signal is deterministic and exactly linear then d2s(x)/dx2 = 0. To allow more flexibility in the model, it will instead be assumed that the signal is stochastic and d2s(x)/dx2 = u(x) where u is a white noise process. So, one can write the following continuous state space model y(x) = (l ,0)a(x) + cre(x), da(x) = Aa(x) + HdWx, and H = d i a g ( T i , T 2 ) . The discrete version of this model is easily seen to be y(xi) = (l ,0)a(x;) + cre(x) a(xi+1) - ehiAa(xi) + v(xi) where hi — — eh'A — Var(v(xi)) — For the cycle model, if one assumes that the signal behaves exactly as a trigonometric function with frequency u>, corresponding to a period p = 2TT/UI, then s satisfies the differential equation s"(x) + UJ2S(X) — 0. Hence, to allow for departures we will assume where A 0 1 0 0 1 hi 0 1 i — 2 , . . . , n, and 76 that s satisfies the stochastic differential equation d2s(x)/dx2 + LO2S(X) = u(x), where u is a white noise process. So, one can write the following continuous state space model where A 0 1 2 0 —LO Now note that if P y(x) = (l,0)a(a:)+ <re(x), da(x) = Aa(x) + H dWx, and H = d iag(r i , T2 ) . ! Q W J , then P-'AP 0 u -LO 0 (3.36) So, one can transform (3.36) by y(x) - Z*a*(x) + ae(x), da*(x) = A*a*(x) + H*dWx, where a*(x) = P~la(x), A* = P ~ M P , H* = P~XH, and Z* = (1,0)P. Note that ehiA* is the same transition matrix as in the discrete cycle (or seasonal) model derived in Chapter 2. The matrix H* is diagonal. Examp le 3.5 (A linear trend plus cycle continuous structural model.) Let us consider the data from Andrews and Herzberg (1985) on the age-adjusted incidences of melanoma for the years 1936-1972. This data set has been analysed, among others, by Heckman and Ramsey (1996) using L-splines. The differential operator used is L = u2D2 + D4. This corresponds to the model of Example 3.3 with T\, T2 and T3 set to zero a priori of the analysis. 77 Heckman and Ramsey compared the L-spline fit with the one given by a polynomial smoothing spline of order 3 (a cubic spline) and 7. A smoothing polynomial spline of order 2m — 1 is an L-spline with differencing operator L = Dm. The comparison of the fits is done in terms of the generalised cross validation criterion as a function of only the smoothing parameter. The fit that gives the lowest generalised cross validation score is chosen as best. In Figure 3.5 one can see the predicted signal given by using the mentioned L-spline and a polynomial smoothing spline of order 7. Figure 3.7 shows the cubic spline fit to the data. The solid lines are the predicted signals and the dashed lines are the predicted error bands. A dot represents a datum. The polynomial smoothing spline is obtained from the following continuous state space model y(x) = Za(x) + ae(x), da(x) = Act(x) dx + r dWx, with Z = (1,0) or Z = (1,0,0,0) for the cubic or polynomial spline of order 7 respectively, e a standard white noise process, W a standard Wiener process, and 0 1 0 0 " 0 0 1 0 0 0 0 1 ' 0 0 0 0 . for the cubic or polynomial spline of order 7 respectively. From this continuous model one can easily derive the following discrete S S M , yi = Za.i + ati, = Ta>i + TUi, 78 A = 0 1 0 0 or A = T and Var(ui) = 1 1 0 1 1 j 2 L or 1 1 1/2 1/3! 0 1 1 1/2 0 0 1 1 0 0 0 1 or Var(ui) 252 72 30 24 X i _ i i 72 20 X I 30 8 1 1 L 24 6 for the cubic or polynomial spline of order 7 respectively. 3 2 The estimated parameters obtained by minimising the generalised cross validation criterion are given in Table 3.2. Model LO a2 T2 G C V L-spline L = LO2D2 + D4 0.650 0.061 0 0.0067 cubic spline L = D2 — 0.107 0 0.0078 Pol. spline order 7 L = D4 0.049 0.014 0.0098 Table 3.2: Estimated parameters for the L-spline and polynomial spline of order 7 models for melanoma data Figure 3.6 shows the predicted innovation residuals (solid lines) together with point-wise error bands (dashed lines) joined by straight lines. Figure 3.7 shows the fit obtained by using a smoothing cubic spline together with con-fidence bands. As one can see the predicted signal is a line. This is judged unsatisfactory since it is known from the reference source of the data set that there is a periodicity in the number of cases of melanoma along years since the disease is associated with the number of sunspots which is a periodic variable. In fact, when plotting the innovation residuals, 79 Figure 3.5: Predicted signal and pointwise confidence bands joined by straight lines of incidence of melanoma data using an L-spline and a polynomial spline of order 7. which if the model is nearly correct should be independent, it appears as if a periodic component is missing in this model. Heckman and Ramsay concluded that the L-spline provides a better fit than the poly-nomial spline of order seven fit based on the fact that the minimum of the generalised cross validation criterion, as a function of the smoothing parameter, is lower for the L-spline fit than for the polynomial spline of order 7 fit. In Table 3.2 the G C V scores were computed but as a function of all parameters involved. As one can see, the G C V score is indeed lower for the fitted L-spline model but the difference between the scores is not big, as was the case in Example 3.4. Moreover, the G C V score for the cubic spline fit is lower than for the polynomial spline of order 7. In Figures 3.5 and 3.6, one can see that the error bands for the predicted signal using 80 L-spline Polynomial spline of order 7 1940 1950 1960 1970 year Hi 1940 1950 1960 1970 year Figure 3.6: Predicted innovation residuals of incidence of melanoma data with error bands for the L-spline and polynomial spline of order 7 models. the polynomial spline model of order 7 are rather wide and the innovation errors are definitely larger than the innovation errors obtained using the L-spline model. Note that the use of diagnostic tools can help the data analyst in assessing the ade-quacy of a fit. The G C V score does not seem to be useful in this case. Also, if one had used the approach by polynomial splines based just on the minimisation of a criterion then one would have not been able to incorporate observed and/or known features of the data, in this case, a trend and a seasonal component. It needs to be stressed that an analysis of the cubic spline model where residuals are available reveals the fact that a seasonal component is missing in the model. In this case we know a priori that the seasonal com-ponent should be present but in many cases one does not know and the methodology for data analysis should allow for the discovery of this kind of feature. One may incorporate these features by means of explicitly expressing the behaviour of the signal as a stochastic 81 Sample Autocorr. Figure 3.7: Predicted signal and pointwise confidence bands joined by straight lines of incidence of melanoma data using a cubic spline and autocorrelation function of the innovation residuals. differential equation or use a structural model (trend + seasonal). • 3.6 Summary In this chapter the concept of continuous state space models has been reviewed. We have dealt with the issue of estimation and prediction using a continuous state space model. For this purpose we have dealt with the derivation of a discrete model starting from a continuous state space model, particularly with the computation of the state error matrix and initialisation of the diffuse Kalmaii filter. To illustrate the increased usefulness of continuous state space models over the non-82 parametric technique that uses L-splines to solve smoothing problems, two examples were analysed. The advantages of a model based approach to smoothing was demonstrated. 83 Chapter 4 Smoothing polynomial splines Researchers working on regression problems realized that global models, such as the ones proposed in parametric regression, failed to perform well on many occasions. Alternative methodologies to parametric regression were developed based on the idea that instead of assuming the signal to be a known, up to some parameters, it can be assumed to be deterministic but of an undetermined form. It is assumed that the signal has some smooth properties such as derivatives of some order, etc. Among these methods is polynomial spline smoothing, where piecewise polynomials come to play an important role. It is claimed that because of their segmented nature, piecewise polynomials are more flexible than (global) polynomials and thereby adjust more effectively to the local characteristics of data. Smoothing polynomial splines are obtained by minimising a penalised least squares criterion where the penalty is the integral of the square of the derivative of some order of the unknown regression function. If the penalty is the integral of the square of the second derivative of the regression function, then the smoothing polynomial spline is called a cubic spline and it is a cubic piecewise polynomial. Cubic splines have attracted much attention as an alternative to parametric regression. 84 Perhaps their popularity is due to the fact that the criterion giving rise to cubic splines can be seen as a combination of a goodness of fit criterion and a term that penalises on excessive Euclidean curvature of the regression function. Also, this criterion can be seen as a penalised least squares criterion where the penalty is a product of some matrices with special characteristics. This enabled algorithms that facilitate the computation of cubic splines. Green and Silverman (1994) dedicate almost an entire book to reviewing cubic splines, computational algorithms and possible uses and applications of this nonparametric methodology. In this chapter a review of the most salient ideas and results in polynomial spline smoothing are presented. At this point, it is worth noting that, as already pointed out in Chapter 3, smoothing polynomial splines are a special case of L-splines. That approach in this chapter is rather different and the following is not based on stochastic differential equations but on the construction of a basis for the space of (natural) splines. This approach is the one originally pursued. In Section 4.1 the definition and a result on the existence of smoothing polynomial splines are given. In Section 4.2 a basis for the vector space of smoothing polynomial splines is presented. This section is mainly devoted to proving the new results stated in Lemma 2 and Theorem 4.2 which will be crucial for the development of the next chapter. These results generalise (to polynomial splines of any possible order) the result that asserts that cubic splines can be obtained as the solution of a certain penalised least squares problem (Green and Yandell, 1985). In Section 4.4.2, a review of work by Wahba (1978) and of Wecker and Ansley (1983) in terms of relating polynomial smoothing splines with stochastic models is reviewed. For a source of material on the topic, see Eubank (1988) and the references therein. 85 4.1 Definition of polynomial smoothing splines For the sake of simplicity assume that the indexing variable takes values on [0,1]. Let L12 [0,1] be the set of all square integrable functions on the interval [0,1] and let the rath order Sobolev space be = {/ : is absolutely continuous, j = 0 , . . . , ra - 1, and / ( m ) £ L2[0,1]} . (4.1) Now, suppose that g £ W^tO, 1]. To estimate g, an estimator is desired which fits the data well but, at the same time, has some degree of smoothness. A natural measure of smoothness associated with a function / £ W 2 m[0,1] is fo[f^m\t)]2dt, while a standard measure of goodness-of-fit to the data is n~lYJj-\(Vi ~ f(xi)Y- Thus, one method of combining these two criterion functions to obtain an overall performance measure is to assess the quality of / by its values for (1 - tin'1 Y.{y3 - f(x3)f + q l\f^\t))2dt i=i J 0 for some 0 < q < 1. A n optimal estimator could then be obtained by minimizing this functional over / £ W ^ O , 1]. However, by setting A = q/(l — q) we see that this is equivalent to estimating g by the function that satisfies the criterion ™ n n _ 1 £ ( w _ ffaW + A f-U^Wdt> A > 0 - ( 4- 2) 2 ' j = l The minimizing function, f\, is the smoothing spline estimator of g. The parameter A in (4.2) governs the tradeoff between smoothness and goodness-of-fit. When A is large (q near 1) this places a premium on smoothness and heavily penalizes estimators with large rath derivatives. Conversely, a small value of A (q near 0) 86 corresponds to more emphasis on goodness-of-fit. These properties have prompted A to be called the smoothing parameter. D e f i n i t i o n 4.1 A spline of order r with knots at x \ , . . . , x n ; Xi < . . . < x n , is a function such that i) it is a piecewise polynomial of order r on any subinterval [ X J , x , + i ) , ii) it has r — 1 continuous derivatives and iii) it has a rth derivative that is a step function with jumps at x i , . . . ,xn. Thus, a spline is a piecewise polynomial whose different polynomial segments have been joined (tied) together at the knots x i , . . .xn. Notice that any piecewise polynomial of order r with r, rather than r — 1, continuous derivatives would actually be a single polynomial of order r. This gives the conclusion that a spline is the smoothest possible piecewise polynomial which still retains a segmented nature. It follows from the previous definition that a spline of order r with knots at x n = ( x i , . . . , x r a), sr>xn, can be written as r n sr,x"{x) = X] ®iX* + &i(x ~ XJY+I £j £ -R, * = 0 , . . . , r, j = 1 , . . . re. i-0 j=l In the previous equation, ( / (x) ) + equals f(x) if / (x) > 0 and zero otherwise. If we denote by <SVjX" the set of all splines of order r with n knots at x 1 ; . . . , x n then it can be easily seen that <SV)X" is a vector space of dimension ra + r + 1, i f x i > 0 and x„ < 1. Of particular importance is the set of natural splines of order r = 2m — 1 with n knots at X \ , . . . , x n . 87 D e f i n i t i o n 4.2 A spline function is a natural spline of order 2m — 1 with knots at xi,...,xn if, in addition to the properties implied by i)-iii) above, it satisfies, iv) it is a polynomial of order m — 1 outside of [xi, xn]. The name natural spline stems from the fact that as a result of iv), the natural spline s satisfies the so called natural boundary conditions 5(J')(0) = s ( j ) ( l ) = 0 , j = m,..., 2m - 1. Let AfSr,xn denote the set of all natural splines of order r with n knots at xi,... ,xn. Then it can be seen that AfSr>xn is a sub-vector space of <SV,X" of dimension n. The importance of natural splines is derived from the fact that a natural spline of degree 2m — 1 with knots at x\,...,xn is the function that minimizes (4.2) over all functions in W ^ f O ; !]• If 2/i5 • • • > J/n cannot be interpolated exactly by some polynomial of degree less than m, then the solution to (4.2) is unique. The following is the main result on the existence of the smoothing polynomial spline. T h e o r e m 4.1 (Lyche and Schumacher, 1973) LetB = {cti,... ,atn} be a basis ofjVS2m-i,xn and suppose that n > m. For fixed 0 < A < oo there is a unique minimizer of (4-2) in W^™[0,1]. Moreover, the minimizer is in j\fS2m-i,xn and hence has the form £3™=i aiai where a = ( a i , . . . , an)' satisfies (X'X + n A * ) a = X'Y, where $ and X are the n xn matrices with (i,j)th entry 88 Jo and respectively. Using Theorem 4.1, the vector of fitted values of Y would be Y = H(X)Y, where H(\) = X(X'X + nX^X'. By analogy with linear regression, the matrix H(\) is called the hat matrix. 4.2 A basis for the space of natural splines In this section a basis for the space of natural splines will be discussed. This basis was developed by Demmler and Reinsch (1975). We focus on this basis because it will make transparent some properties of polynomial smoothing splines. Before going to the basis itself, it will be necessary to introduce some notation on divided differences. 4.2.1 Divided differences Let xi,... ,xn be real numbers such that a;,- < Xj, if i < j. 89 D e f i n i t i o n 4.3 The mth order divided differences of a function g are defined recursively by DL(g) = " P ° " ' ^ ~ g " * ^ . « = m + 1 , , , , , n ; •^ m-t-t—1 1 uwift (^flo = flf(»i), * = l , . . . , n . Tfte m i / i order normalized divided differences of a function g is NDlm(g) = (xm+i-i — Let aJ-)Trl, j = 0 , . . . , ra, denote the coefficient of g(xi-j) in ND'm(g). E x a m p l e 4.1 (First and second order divided differences.) Let hi = Then the first divided difference of g is • g(xj) - g(xi.1) ^1(0) = ^ > t = 2 , . . . , n . The second order divided difference is £>2(</) = (hi + hi+i)-1 If ra = 1, a 0 1 = 1/hi and ct^ = —1/hi. If ra = 2, a 0 2 = 1/ft;, = — (I/hi + l / f t , _ i ) , and a 2 2 = l /ft , -_i . Note that if the a;,-'s are equally spaced with ftj = ft, i = 2 , . . . , re, then c/. = «• = - ( mWv' « j , m - aJ,m h \ j J • 90 ft,-^ ^ (ft; ft — J fl'(si-i) + T— 9(xi-2) , i = 3, , n . 4.2.2 The Demmler-Reinsch basis The basis for natural splines that will be presented has been studied by Demmler and Reinsch (1975). Let Qm : R2 —• R be the function defined by Qm{s,t)= j\s-u)™{t-u)™du, 0 < M < 1 , (4.3) and the associated n xn matrix Qm with (i,j)th entry Qm(xi,Xj). We will see later that Qm is the covariance kernel of a certain Gaussian process called m-fold integrated Wiener process 1 . As a consequence, the matrix Qm is positive definite. Also, let T be the n x m matrix with (i,j)th entry [T],j = a ^ - 1 , i — 1 , . . . , re, j — 1 , . . . , m and A T O be the (n — m) x n matrix, such that the i th row of A m is ( A m ) 4 = ( 0 , . . . , 0 , < ™ , . . . , < ™ , 0 , . . ^ ) ' , i = l,...,n-m. (4.4) i—1 n—m—i Note that AmT = 0. Let ftm = AmQm-iA'm, (4.5) and Fm = ( T ' Q T O - I ^ 1 ) - 1 ^ 1 ^ ™ - ! - Demmler and Reinsch define the basis for J\fS2m-i,xn to be B — {OJI, . . . , an} with (cvi(x) , . . . , a n (x)) = (1, x,..., x m _ 1 , Q m _ i ( x , x j ) , . . . , Q m _ i ( x , xn)) where E is the nxn matrix whose columns are the orthonormal eigenvectors of A ^ O " 1 A m . The first m columns of E correspond to the eigenvectors of the m zero eigenvalues of ^ h i s process will play a central role in defining a Bayesian stochastic model for spline smoothing. 91 FmE A'mn^AmE A^fip/Am and the other n — m columns correspond to the eigenvectors of the positive eigenvalues, in ascending order. The choice of this basis for A/c>2m-i,x" allows to express the mth derivative of a natural spline in a convenient way, which makes possible to translate the problem of finding the smoothing spline estimator of a function to solving a penalised least squares problem. L e m m a 2 If f £ MS 2 m - i , x n then fQ[f{m){u)]2 du = a'A'm^Ama, where a = ( / ( x x ) , . . . , f{xn))'. F •*• m A' n^A. mm T (4.6) P r o o f . One can show that / (x) - (1, X, . . . , X m _ 1 , Qm-l(x, X i ) , • • • , Qm-l(x, Xn)) Also, (Qm-i(x,xo,...,cm_x(x,xB)) = vixyr + p~Tiy ((x-Xi)+m_1'• • • ' ( x - ^ ) l m _ 1 ) . (4.7) where v(x) is a m x 1 vector that involves the x's but not the x,'s. From (4.6) and (4.7) and the fact that AmT = 0 we have that / ( m ) ( x ) = ~^~Tj\ ({x ~ • • •'{x ~ A-°™Ama' (4-8) so that i\fm\u)fdu = o , A / m n - 1 A m g m _ 1 A , m n - 1 A m a (4.9) Jo = a 'A' m f2- 1 A m a. (4.10) 92 • The following result shows how the well known result that cubic splines can be ob-tained by minimising a penalised least squares criterion (see e.g. Green and Yandell, 1985), where the penalty is the product of some special matrices, can be generalised to polynomial splines of all orders. This result, as will be shown in the next chapter, will make transparent some properties of polynomial splines. Theorem 4.2 Criterion (4-2) is equivalent to mm{Y - a)'(Y - a) + n\a'A'Jl^ Ama. (4.11) Proof. Obvious from the preceding Lemma. • Next, some special properties of the matrix S l m in (4.11) will be shown. Definition 4.4 A symmetric n x n matrix A = [a t j ] j ' j=i,.. . , n is said to be 2k + I banded if a>ij = 0 whenever \i — j\ > k. Lemma 3 The matrix ttm is 2m — 1 banded. Proof. The (i, j ) th entry of S l m is m m r=0 ;=o / i \2 n [m i r m • I 1 / V ^ „ « + m / ' ~ „ . \m—1 „ J + " i / „ „. \m—1 = [jm-^iy.) Jo IS r 'm i l + r ~ ) + J[ga'.« ^+'-")+ = ^ m ]_ xy^j JQ Ni,m(u)Njtm(u)du, 93 where m r=0 = NDkm{{x-u)^). Note: • Nk,m(u) = 0 if U > Xk+m Or U < Xk-This follows from the fact that (a; — is a polynomial of order m — 1 in x if x > u and that the mth divided difference of a polynomial of order m — 1 vanishes. • Ni>m(u)Njim(u) = 0 if \i - j\ > m - 1. From the previous remark, Nitm(u)Njtm(u) ^ 0 if and only if u £ [a:;, x , + m ) U [xj,Xj+m). But since [xi,Xi+m) U [xj,Xj+m) = 0 if |i — j\ > m — 1, the remark follows. The proposition follows easily from the last two remarks. • As an immediate consequence of the remarks in the proof of the previous lemma one can easily see that if the x's are equally spaced then the diagonals of the matrix 0 m are constant. E x a m p l e 4.2 (The matrices O i and Q,2) 94 O i = I. Also, 0 2 = (h2 + h3)/3 h3/6 0 h3/6 ••: ••. 0 ••• '•• 0 0 K/6 0 hn/Q (hn-i + hn)/3 (4.12) If hi = h, i — 1 , . . . , n then VLo = h 2/3 1/6 0 1/6 •-. ••. o ••• '•• 0 ; ••. ••. ••. i / 6 0 . . . 0 1/6 2/3 • 4.3 Estimation of the smoothing parameter Cross-validation seems to be the most popular method to choose the smoothing parameter. The basic principle of cross-validation is to leave the data points out one at a time and to choose that value of A under which the missing data points are best predicted by the remainder of the data. Craven and Wahba (1979) propose the use of a related method called generalized cross-validation. In this thesis, a review of the cross validation estimation method was presented in Chapter 2. 95 4.4 A Bayesian approach to spline smoothing In this section, a pioneer idea by Wahba (1978) and a subsequent of Wecker and Ansley (1983) that relate polynomial splines with stochastic processes are presented. Wahba shows that polynomial splines can be obtained by setting a signal plus noise model where the signal has a prior distribution which is the same as the distribution of a special kind of Gaussian process, an integrated Wiener processes of some order. Then she shows that the polynomial smoothing spline is the conditional expectation of the signal given the variable that is observed (posterior mean of the signal). Wahba's model has the disadvantage that it does not suggest an algorithm for computation. Later, Wecker and Ansley (1983), based on Wahba's model, develop a state space model for smoothing splines of all orders. Although Wecker and Ansley's model suggests the use of Kalman filtering and related tools for estimation and prediction, the fact that Wahba's model has diffuse initial conditions posed a problem which is dealt with in quite a cumbersome manner. This section has a merely informative purpose. In particular, it will be shown how the matrix f i m in criterion (5.3) can be seen as the covariance matrix of a particular Gaussian process. 4.4.1 Some Gaussian processes Def in i t ion 4.5 Let us denote a stochastic process in [0,1] by {Y(t) : t € [0,1]}. The process is said to be Gaussian if {Y{t\),... ,Y(tr)} has a multivariate normal distribution for any finite collection of points ti,...,tr from [0,1]. 96 To specify a Gaussian process it suffices to specify its mean value function u(t) = E[Y(t)} and its covariance function or covariance kernel K(s,v) = Cov (Y(s),Y(v)). Examp le 4.3 (White noise processes.) A white noise process is a Gaussian process with mean zero and covariance kernel Examp le 4.4 (Wiener processes.) Another example of a Gaussian process is the Wiener process W = {W(t) : t € [0,1]}, which satisfies • W has independent increments, that is, if u < v < s < t then W(t) — W(s) and W(v) — W(u) are independent; • W(s +1) — W(s) is A/"(0,cr2r) for all s,t > 0 where a2 is a positive constant. The covariance kernel of W is given by c(s,t) = a2 min{s,t}. W is called a standard Wiener process if a2 = 1. (4.13) The white noise is said to be standard if a = 1. • 97 • Examp le 4.5 (Integrated Wiener processes.) Let W denote a standard Wiener process and let Wm(t)= / * ^—dW{u), 0<t<l, m = 0 ,1 ,2 , . . . (4.14) Jo m! denote the m-fold integrated Wiener process. Then, W 0 (t) = W(t), Wm(t) = Wm-i(u) du, m — 1,2, The process Wm has covariance kernel Qm, as defined in (4.3). • Proper t ies of integrated W iene r processes Next some properties of integrated Wiener processes will be shown. These properties will be useful to establish several stochastic models to obtain smoothing splines as well as to determine the relations and equivalences of these models. A property of m-fold integrated Wiener processes will now be shown. It is used in Wecker and Ansley (1983) to develop a state space model for spline smoothing based on Wahba's Bayesian model (see Section 4.4.2). In particular the following result gives a dynamic linear representation for m-fold Wiener processes. L e m m a 4 Let wm(x,y)= I" {V ~ ^ dW{u). Jx ml Then, for x,y,z 6 [0,1], m-1 (y — zy wm(x, y) = wm(z, y)+^2 H wm-j(x, z). (4.15) j=o J • 98 Proof. This proof is taken from Wecker and Ansley, 1983. Jx m l Jz m l = wm{z,y) + ^ J—\ J-*-dW{u) Jx ml ™ (y - z)j f* (z - U)m~j ,„., , i=o J • • The following property characterizes the properties of an integrated Wiener process when observed on a discrete grid. Lemma 5 Let Wm denote an m-fold integrated Wiener process. Then the variance matrix of Wm is Qm, where Qm is given by (4-3). Moreover, the covariance matrix of the mth order divided difference ofWm-\ over a grid x i , . . . , x n is f i m - i , as defined in (4-5) and is hence 2m — 1 banded. Proof. That the variance of Wm is Qm follows from Theorem A.1 . The other assertions follow from the definition of fim and Lemma 3. • 4.4.2 Wahba 's approach Suppose that the observations ( x i , yi),..., (xn, yn) are obtained by sampling from a stochastic process [Y(x) : x £ [0,1]} at predetermined points X i , . . . , x n . The process Y 99 is assumed to have a special form with Y(x) = n(x) + t(x), x £ [0,1], (4.16) where {e(x) : x £ [0,1]} is a white noise process with var(e(x)) — a2. A Bayesian structure is introduced by placing a prior distribution on the mean value function fi of the process Y. It is assumed that {fi(x) : x £ [0,1]} has the same prior distribution as the stochastic process m-l j J2 0j- + rZ(x), r > 0 , (4.17) j=o 3 • where 0 = (Po, • • • , P m - \ ) ' n a s a zero mean, normal distribution with Var(P) — vl and {Z(x) : x G [0,1]} is an m - 1-fold integrated Wiener process. The random variable P and the process Z are assumed to be uncorrelated. In model (4.16)—(4.17), instead of treating fi as a deterministic function in the mth order Sobolev space, it is treated as being a polynomial plus a special type of random error function. A possible interpretation of the Bayesian model given by (4.16)-(4.17) is that we extract a polynomial trend from /x and then approximate the remainder by a stochastic process. In fact, the actual role of the Pj is to allow Wo to start at Po instead of at zero. Wahba (1978), showed that the smoothing spline rin,\{-) which is the minimizer of (4.2) has the property ^ x ( x ) = \lmEu{fi(x)\Y = y} (4.18) with A = a2/TIT2, where Eu denotes the expectation under the posterior distribution of /j,(x) with the prior (4.17). (u — oo corresponds to a "diffuse" prior on 0). One important theoretical drawback of model (4.16)-(4.17) is that it is not equivalent 100 to the model Vi = Kxj) + e j , 3 = 1, • • • ,n, (4-19) where a < x\ < . . . < xn < 6, e i , . . . , e n are mutually independent and identically distributed with mean 0 and variance a2 and fx G W™[0,1]. If // in (4.16) and (4.17) had an rath derivative it would have to be the derivative of a Wiener process which is known to be a smooth process that is nowhere differentiable with probability one. Hence, the sample paths of the process (4.17) have null probability of being in the rath order Sobolev space. However, Ltn,\ defined in (4.18) is most definitely a member of this function space (see Wahba, 1983). Consequently, the Bayesian approach gives an estimator of LI which is not a plausible candidate for the realized regression curve. This estimation paradox causes some conceptual difficulties. 4.4.3 A state space form of Wahba's model This is the approach to spline smoothing developed by Wecker and Ansley (1983). The authors based their work on the stochastic model for LI(X) given by Wahba (1978), as in (4.17) with a diffuse prior on 0. A state space model is derived to estimate the parameters involved in the model as well as to construct an estimate for fi(x) together with confidence bounds. The authors make an extensive use of property (4.15) to construct a state space model. If in (4.15) we let x = X\, z = x ; _ i and y = we have that m hj-Wm(xuXi) = U > m ( 3 i - l , 3 j ) + 2 ^ l W m _ j ( & i , X i _ i ) , (4.20) j=o 3-with hi = X i and hi = Xi — x 8 _ i , i = 2 , . . . , n. Define OLi = T ( w m - i ( x U X i ) , W m - 2 { x U X i ) , . . . , W o ( x U X i ) ) ' 101 Ui = T ( w m - 1 ( x i - i , X i ) , W m - 2 ( X i - l , X i ) , . . . , W o ( X i - 1 , X i y . Let T(u) be the m x m matrix T(u) 1 u 0 1 0 0 ( m - l j ! u m - 2 (ro-2)! It is easy to see that T(u + v) = T(u)T(v). Then, (4.20) implies that O ; = T,OJ;_I + I f ; , where T; = T(hi). and the variance matrix of u;, 0;, has jkih entry [Sli]jk = r (xi-h)m-i(xj-h)m (m — j)\(m — k)\ ^2m-(j+fc)+l dh — T (2m - (j + k) + l ) (m - j ) ! (m - k)\' (4.21) Letting / = ( 1 , 0 , . . . , 0)', ft — (fto, • • •, ftm-i)', Wahba's model can be written in state space form: yi = Xift + I'cti + €.{ ai = T ; a ; _ i + Uj , i = l , . . . , n , (4.22) where X{ = (1 ,£ ; ,xf/2\, . . . ,x™ 1 / ( m — 1)!) and starting conditions ao = (0 , . . . ,0) ' . It follows from the definitions that a ; _ i is independent of u; . The variance of u ; is given by f i ; . Note: model (4.22) implies that not only the state vector and its variance parameters but also the vector of parameters (fto, • • • ,ftm-i) should be estimated. 102 At a first glance (4.22) looks like a semiparametric model for ( ? / i , . . . , yn). This is misleading because the matrix X = (X[,..., X'n) cannot be interpreted as a covariate and the parameters 0 do not have an obvious interpretation as they stand in the model. The reason why these parameters and this matrix appear in the model is the following: an m — 1 integrated Wiener process has the constraint that its value at zero as well as its m — 1 derivatives are the zero random variable. B y adding the term X0 to the state one allows for nonzero initial conditions for the process and its derivatives (the sum of X0 and Wm is known as a free m-fold Wiener process because it allows for free initial conditions). Now, observing that Xi = l'T(xi) one can rewrite (4.22) as Vi = l'(T(xi)0 + ai) + ei T(xi)0 + ai = T(hi + x,-_i)/3 + r ( / i t )a i _ 1 + m = T(hi){T(xi-1)0 + at-_x) + u{. If we define a* = T(xi)0 + ai we observe that y{ = l'a* + 6i , i = l,...,n, (4.23) a* = Tia*_x + Ui with starting conditions Q;Q = 0. Therefore, if Wahba's model is to be stated in state space form, the formulation given by (4.23) should be preferred. Note that the D K F provides a transparent and easy to use prediction and estimation tool for this situation. 103 4.5 Summary In this chapter, the concept and some properties of polynomial smoothing splines were reviewed. As it has been shown in Chapter 3, polynomial smoothing splines are a special case of L-splines. Nevertheless, the approach in this chapter is the originally pursued. A basis for the vector space of smoothing polynomial splines was reviewed with the purpose of making some properties of polynomial smoothing splines more transparent. In particular it is shown that the penalised least squares criterion that polynomial smoothing splines minimise is equivalent to another such criterion, with a penalty given by the product of some matrices. This is crucial for the development of the next chapter. For the sake of completeness of the review of polynomial splines, some work that connects polynomial splines with stochastic processes is presented. 104 Chapter 5 A closer look at polynomial splines Smoothing by cubic splines has become in the latest years a popular alternative to para-metric regression. The popularity of this smoothing methodology is due perhaps to the penalised least squares criterion giving rise to it. This criterion is a linear combination of an ordinary least squares criterion and a measure of global roughness or total amount of Euclidean curvature of the regression curve. A regression curve minimising such a crite-rion should provide both fidelity to the data and a good degree of smoothness. Indeed this combination of criteria is appealing and might seem infallible for regression data analysis. Nevertheless, as it has been shown in Chapter 3, Example 3.5, a cubic spline does not always provide an accurate or satisfactory description of the data. The scientific mind will immediately question why this is happening. In order to answer such an inquiry, focusing only on the criterion giving rise to cubic splines one might not find an answer. Also, in light of Example 3.5, one may question if it is possible to know more generally, and not just for a particular case, when cubic splines might be an appropriate regression methodology to use or not. As will be shown, a cubic spline can be seen as an A R I M A ( 0 , 2 , 1 ) plus noise process and so, for example, periodic data will unlikely be well described by cubic 105 splines. This conclusion is arrived at by looking at a stochastic model that gives a cubic spline; searching for information about cubic splines beyond the penalising criterion that gives rise to them. In this chapter, based on the material presented in Chapter 4 and some new results, smoothing polynomial splines will be explored with more detail. First, it will be proved that the penalising criterion that defines the smoothing polynomial splines is equivalent to a penalised least squares criterion where the penalty is a quadratic form. B y stating the solution of a penalised least squares problem into a signal extraction problem one can state the solution to the cubic splines problem as a signal extraction problem as well. Next it is shown that the implications of the previous connections is that a model behind polynomial spline smoothing of order 2m — 1 is an A R I M A ( 0 , m, m — 1) plus noise, with all coefficients fixed, except for the noise variance. 5.1 T h e s igna l p lus noise m o d e l Suppose that observations (yi,X{), i = l,...,n, are obtained. Further, assume that the observations consist of a signal (the variable of interest) plus a noise. Let Y denote all the observations made, s the unobserved signal and e the noise. The signal plus noise model is the relationship Y = s + e, (5.1) where E(s) = E(e) = 0, Var(e) = XI, Cov(e,s) = 0 and Var(Y) = E . If the signal plus noise model holds we have that Var(s) = E — XI and Cov(s,Y) = Var(s). The best linear predictor of s given Y, P(s\Y), as shown in Chapter 2, is P(s\Y) = Cou(s,r)[Var(y)]- 1 r 106 = Varis^Var^ + Varte)}-^ = [7 + y a r ( e ) ( y a r ( 5 ) ) - 1 l " 1 r . Then, it follows that the predictor of the signal, s, is s = {l + \[Var(s)]-iy1Y. (5.2) For a proof of 5.2 see, e.g., Green and Silverman (1994), pp. 18-19. The estimator s is called the statistical smoother. 5.2 Penalised least squares Let Y and s be n-dimensional vectors and A be a n x n , full rank matrix, called the penalising matrix. The penalised least squares problem consists of finding s G Rn that satisfies m m ( y - s)'(Y - s ) + Xs'As. (5.3) The solution to the penalized least squares problem is 5 = (7 + A 4 ) _ 1 F . (5.4) Observe that proper choices of A makes penalized least squares suitable for smoothing. As we shall see, the penalized least squares criterion used to find the smoothing spline (of any order) can be put into a penalized least squares criterion like (5.3). Hence, to translate the penalized least squares problem into a state space statistical formulation, one can construct a model such that Y = s + e and Var(s) — A - 1 and perform a signal extraction in this context. 107 5.3 Smoothing splines, penalized least squares, sig-nal extraction and state space models From (5.2) and (5.4), the penalized least squares solution is the statistical smoother if and only if Var(s) = A - 1 . Theorem 4.2 implies that spline smoothing can be viewed as a penalized least squares problem wi th penalising matr ix A ^ f i " 1 Am. L e m m a 6 If Var(Ams) = flm then [Var(s)}-1 = A'mil^Am. P r o o f . Let ( A m ) £ and ( A ^ ) ^ denote a left and right inverse of Am and A'm respectively, which exist since rank of A m = n — 2. If Var(Ams) = flm then fl. = AmVar{s)A'm. Then, Var(s)A'mCl~1A. ( A m ) l Z A m I / a r ( s ) A , m O - 1 A ( A m ) £ O m O " 1 A and A ^ 0 " 1 A m y a r ( 5 ) A' (T 1 A Var(s)A'm(A'J-R A ^ f 2 - 1 0 m ( A ' m ) ^ I. Hence, [Var(s)}-1 = A ' ^ A . • 108 Based on the preceding Lemma, the next result shows how a signal plus noise model can be constructed, with a stochastic signal, so that the solution to the problem of signal extraction or prediction is the same as the polynomial smoothing spline of order 2m — 1. T h e o r e m 5.1 If Y = s + e, and Var(Ams) = Clm, then the solution to the signal extraction problem is also the solution to the problem of finding the polynomial smoothing spline of order 2m — 1. As will be seen, Theorem 5.1 can be used to formulate a state space model, bringing not only computational advantages but also making more transparent the properties of smoothing splines. Following we show how to construct a state space model for smoothing with cubic splines. E x a m p l e 5.1 (A state space model for cubic spline smoothing) Let r / i , . . . , yn and a i , . . . , an be random variables such that (5.5) where i — (^8 5 ^ *t) ) ^ = l , . . . , n ; 7V(0,cr 2); 109 • u i , . . . ,un are independent random variables and u,- ~ N2(0, E,) where / U ? + |/z 2 + fefc,- \ ^ ~ ^ \h} + ^ Zl,- J where k{ is any number. • ui , . . . , u n are independent of OJI, . . . , a n . Next it is shown that, according to the model specified in (5.5), Var ( A 2 s ) where H 2 is as in (4.12). Let T i = { 0 l ' ) ' 1 = 1 ' - - - ' n -From (5.5) we have that ai — ai-i = (Ti — I)ai-i + u,-. Also, note that if we define 1 (1 1 \ 1 oti = —ai — I - — h T a,-_i + a,-_2 = (ai - ai_i) /hi - (at-_i - a,-_2) = ft-^Ii - J)ai_i - hj2x(Ti - 7)at_2 + ft~V- - /V-i^-i then A 2 a = ( d 3 , . . . , an)'. Letting s,- denote the first component of d,- and u; = (u|,u 2 ) ' we have that Si = rt-_i - r;_2 + / i t _ 1 u- - / V J W - I -Also from (5.5), r; = r t _i + u 2 = r 0 + u 2 + • • • + u2. Therefore, Si = w2_x - V-Vi- i + W -110 Simple calculations give Var(si) = (hi-i + hi) /3, i = 3 , . . . , n Cov(5iiSi+i) = hi-i/6, i = 3 , . . . ,n — 1, Cov(si,si+j) = 0, i f j > 2 , as asserted above. • 5.4 Smoothing splines and A R I M A models In this section it is shown how smoothing with splines is equivalent to modeling the response variable y as an A R I M A model plus noise with fixed coefficients. Let us assume that m — 2 (i.e we want to smooth by means of cubic splines). Let the Xj's be equally spaced with h = x,- — x ; _ i , i — 2, . . . , n . Without loss of generality, take h = 1. Define Si — 2s,_i — S j _ 2 + Ui + dUi-i where u,'s are i . i .d A/"(0,r 2) and T 2 , 0 are such that T2(l + [32) 2/3 (5.7) T20 = 1/6. (5.8) Possible values for 0 are 2 ± \/3 and r 2 = 1/(6/3). Then, st- is an A R I M A ( 0 , 2 , 1 ) process and with the specified coefficients Var(A2s) = 0 2 (see Example 4.2). I l l 1 Ui Now consider the following state space model yi = ( 1 0 )at- + f-i where e,-s are iid j\f(0,a2), e,-s and UiS are uncorrelated, and a,-conditions OJQ = (a^a 2 , ) ' . It is easy to verify that aj = Si. (5.9) (a],a2)' with starting The solution to the signal extraction problem is the solution to the smoothing by cubic splines problem. Now assume that the not equally spaced. Suppose that the z's are recorded up to the Z'th decimal (e.g. 0.053 is recorded up to the third decimal) and let h = 10~'. Construct a new equally spaced grid Xj. For this new grid we will have observed j/j only if XJ = Xi for some i = 1 , . . . , n. Model (5.9) still holds but we will express in terms of ctjl_1 where we have observed and y;. Let P be the matrix P = 1 0 - 1 1 Then, P " 1 = 1 0 1 1 and P - l ' 2 1 " " 1 1 ' P = - 1 0 0 1 T. to Hence, letting a* — P~la and H = (1,0 + 1)' one has that model (5.9) is equivalent V i = (1 OK + £ i OL* — T a * _ ! + Hui. Therefore, a* = T h a t h + J2Th-jHui-h+3, i=i 112 where 1 h - j 0 1 , j = o,... h. Let u. * = zUTh~JHu^- Then> h Var(u*) = r 2Th~JHH'(Th~^)' (5.10) where is a constant. Since (5.10) is of the form given in (5.6) it follows that even in the case when the x's are unequally spaced, the underlying process is an ARIMA(0,2 ,1) plus noise with all the A R I M A parameters fixed. One can repeat the process above if one is dealing with higher order splines and conclude that smoothing with splines of order 2m — 1 is equivalent to modeling the response variable y as an A R I M A ( 0 , m , m — 1) plus noise with all the M A coefficients fixed. Next we show an equivalent model to the one given in (5.5). This equivalent model is more suitable for estimation. Let Pi be the matrix Pi = 1 0 l/k 1/hi Then 1 0 - 1 k and PrlTiPi = 2 1 - 1 0 113 so that (5.5) can be written as y i = ( i O K + ^ 2 1 -1 0 where a* — P a ; a n d u* = P Uj . Hence, F a r « ) = ^ 1/3 + h/hi 1/6 1/6 1/3-ki/hi (5.11) (5.12) 5.5 Related topics 5.5.1 Smoothing in the presence of heteroskedasticity Let H ^ b e a matrix and consider a signal extraction problem in the following model Y = s + We (5.13) where Y, s and e are as in Section 5.1. It is easy to see that the statistical smoother is given by s = {WW + ^ ^ ( j j j - 1 } - 1 ^ . (5.14) Note that s is also the solution to the penalised least squares minimization problem min(r - s)'WW'(Y - s) + Xs'As, s6Rn (5.15) where A = Var(s) - l It is easy to build an appropriate state space model for smoothing in the presence of heteroskedasticity. Simply choose an appropriate matrix G in Definition 2.1. 114 5.5.2 Repeated measurements Suppose that at some point x,- one observed several values of y, y,j, j = 1 , . . . , m,-. Further assume that E(yij) does not depend on j but only on i (hence, the repeated measurements are taken with the same mean). A reasonable penalised criterion that one could use for smoothing is n mi . E E f e r W t A r ' f . (5.16) i=l j=l It is easy to verify that one can write (5.16) as (Y - s)'M(Y - s ) + Xs'Ks; (5.17) with K chosen appropriately, M a diagonal matrix with diagonal entries ra,- and Y rep-resenting the vector with entries y,- = ra"1 YJf=\ yij-5.5.3 Semiparametric regression A useful generalization of (1.1) is the following V i = Xi/3 + g(Xi) + e,-. (5.18) There are a number of situations where this type of model might be appropriate. Note-worthy is the problem of analysis of covariance where P represents a vector of treatment effects and x is a vector of covariates. The classical method of analysis would assume that the covariates are also linerly related to Y. But there are cases where non-linearity in x needs to be contemplated. 115 Looking at the state space model defined in Definition 2.1 we see how naturally one could introduce a covariate X without having to device any new computational methods. So one can easily solve semiparametric regression problems. A n d we can make our choice of fixed or random effects by specifying C appropriately in Definition 2.1. 5.5.4 Change point problems or intervention analysis These are the types of problems where the signal presents one or more discontinuity and which will not be removed by simple differencing. The main idea is to study the effect of certain events or "shocks" that produced a sudden change in the signal. Mueller (1992) analyses the problem from a non-parametric point of view. The result is a quite cumbersome approach to intervention analysis. Harvey (1989, Chapter 7) treats the topic in an illustrative, accessible and interesting way. 5.6 Summary In this chapter, and based on the results of Chapter 4, the connection between polynomial splines and A R I M A processes was discussed. This connection should serve to understand the nature of polynomial spline smoothers. In particular it is shown how this methodology, that is presented in the literature as a powerful tool for regression, is in fact quite limited. 116 Epilogue In this thesis the increased usefulness of the utilisation of state space models as a smoothing tool has been shown. A model based approach to smoothing permits the assessment of the predictions and estimates obtained in a natural, clear and direct way. For each nonparametric technique of smoothing based on the minimisation of a pe-nalised least squares criterion, there is an associated state space model. It turns out that this state space model usually does not have all the flexibility one could have because the focus is put into minimising the penalised least squares criterion rather than on exploring the model at hand. It has been shown that, starting from the same idea that the signal will be close in form to the solution of a differential equation, the continuous state space approach to smoothing results in insightful data analyses. The advantage of posing a smoothing problem as a signal extraction problem in the context of a state space model is that existing efficient algorithms and diagnostic proce-dures become readily available for estimation and model assessment purposes. Also, the researcher is provided with tools to understand the results of the data analysis procedure and to look for possible ways of improving the analysis. In order to compare L-splines with state space models, some important computational results and details about implementation of the prediction and estimation tools had to be derived. As a way of making transparent certain features of polynomial splines, we found that the criterion giving rise to polynomial smoothing splines is equivalent to a penalised least squares criterion where the penalty is a product of matrices. This enables us to show the connection between polynomial splines and A R I M A processes. 117 References Anderson, B . D . O . , and Moore, J .B . (1979). Optimal Filtering. Prentice-Hall, Englewood Cliffs. Ansley, C . , and Kohn, R. (1986). O n the equivalence of two stochastic approaches to spline smoothing. Journal of Applied Probability, 23 A , 391-405. Ansley, C . , Kohn, R., and Wong C . (1993). Nonparametric spline regression with prior information. Biometrika. 80, 1, 75-88. Arnold, L . (1974). Stochastic Differential Equations: Theory and Applications. John Wiley & Sons, Inc. New York. Bellman, R. (1960). Introduction to Matrix Analysis. McGraw H i l l , New York. Bergstrom, A . R . (ed.) (1976). Statistical Inference in Continuous Time Economic Models. North Holland, Amsterdam. Bergstrom, A . R . (1984). Continuous time stochastic models and issues of aggregation over time. In Z. Griliches and M . Intriligator (eds.) Handbook of Econometrics, vol. 2, 1145-1212. Box, G . (1979). Robustness in the strategy of scientific model building, in Robustness in Statistics, edited by Launer, R . L . and Wilkinson G . N . . Academic Press, New York. Box, G . (1980). Sampling A n d Bayes inference in scientific modeling and robustness (with discussion). Journal of the Royal Statistical Society, Series A. 143, 383-430. Box, G . , Jenkins, G . , and Reinsel G . (1994). Time Series Analysis. Forecasting and 118 Control. Prentice Hall , Englewoods Cliffs, New Jersey. Cook, R., Weisberg, S. (1982). Residuals and Influence in Regression. Monographs on Statistics and Applied Probability. Chapman and Hall , New York, London. Crafts, N . F . R . , Leybourne, S.J., and Mills, T . C . (1989). Trends and cycles in British industrial production, 1700-1913. Journal of the Royal Statistical Society, Series A , 152, 43-60. Craven, P. and Wahba, G . (1979). Smoothing noisy data with spline functions. Nu-merische Mathematik, 31, 377-403. de Jong, P. (1987). Linear Time Series for statistical analysis and prediction. Unpublished manuscript. de Jong, P. (1988). The likelihood for a state space model. Biometrika, 75, 1, 165-169. de Jong, P. (1988). A cross-validation filter for time series models. Biometrika, 75, 3, 594-600. de Jong, P. (1989). Smoothing and interpolation with the state space model. J. Amer. Statist. Assoc., 84, 1085-1088. de Jong, P. (1991). The diffuse Kalman Filter. The Annals of Statistics, 19, 1073-1083. de Jong, P., and Chu-Chun-Lin , S. (1992). Stationary and non-stationary state space models. Journal of Time Series Analysis, 15, 151-166. de Jong, P., and Chu-Chun-Lin , S. (1994). Fast likelihood evaluation and prediction for nonstationary state space models. Biometrika. 81, 133-142. 119 Demmler, A . , Reinsch C . (1975). Oscillation matrices with spline smoothing. Numerische Mathematik, 24, 375-382. Downing, D . J . , Pike, D . H . , and Morrison, G . W . (1980). Application of the Kalman filter to inventory control. Technometrics, 22, 17-22. Draper, N . , and Smith, H . (1981). Applied Regression Analysis. Second edition. Wiley, New York. Eubank, R . L . (1988). Spline Smoothing and Nonparametric Regression. Marcel Dekker, inc., New York - Basel. Gallant, A . R . (1977). Testing a nonlinear regression specification: a nonregular case. J. Amer. Statist. Assoc., 72, 523-530. Gardner, G . , Harvey, A . C . , and Phillips, G . D . A . (1980). Algorithm AS154. A n algorithm for exact maximum likelihood estimation of autoregressive-moving average models by means of Kalman filtering. Applied Statistics, 29, 987-1007. Golub, G . H . and Van Loan C . F . (1983). Matrix Computations. John Hopkins University Press, Baltimore. Green, P.J. , Silverman B . W . (1994). Nonparametric Regression and Generalized Linear Models. Chapman & Hall , London. Green, P.J. , Yandell, B . (1985). Semi-parametric generalized linear models. Generalized Linear Models. Lecture Notes in Statistics. 32, 44-55. Springer, Berlin. Harvey, A . C . and Stock J . H . (1986). Estimation, smoothing, interpolation, and distribu-tion for structural time-series models in continuous time. Economic and Social Research Council and The London School of Economics and Political Science, London. 120 Harvey, A . C . (1989). Forecasting, Structural Time Series Models and the Kalman Filter. Cambridge University Press. Cambridge. Harvey, A . C . (1993). Time Series Models, Second Edition. Harvester Wheatsheaf, New York. Harvey, A . C , and Durbin, J . (1986). The effects of seat belt legislation on British roar casualties: a case study in structural time series modelling. Journal of the Royal Statistical Society, Series A , 149, 187-227. Harvey, A . C , and Stock, J . H . (1986). Estimation, smoothing, interpolation and distribu-tion for structural time series models in discrete and continuous time. The London School of Economics and Political S c i e n c e / D E M E I C Econometrics Programme Discussion Paper A62. Heckman, N . E . and Ramsay, J . (1996). Beyond polynomial splines: customizing the spline basis. Paper presented at the meeting "Program on spline functions and the theory of wavelets", "Splines and wavelets in statistics" section, A p r i l 9-12, 1996. Kimeldorf, G . , and Wahba, G . (1971). Some results on Tchebycheffian spline functions. Journal of Mathematical Analysis and Applications. 33, 82-95. Kohn, R., and Ansley, C . (1987). A new algorithm for spline smoothing based on smooth-ing a stochastic process. SI AM Journal on Scientific and Statistical Computing. 8, 1, 33-48. Lyche, T . and Schumaker, L . (1973). Computation of smoothing and interpolating natural splines via local bases. SIAM Journal of Numerical Analysis, 10, 1027-1038. Mueller, H . G . (1992). Change points in nonparametric regression analysis. The Annals of Statistics, 20, 737-761. 121 Phadke, M . S . (1981). Quality audit using adaptive Kalman filtering. ASQC Quality Congress Transactions - San Francisco, 1045-1052. Ramsay, J . (1996). Principal differential analysis: data reduction by differential operators. Journal of the Royal Statistical Society, Series B. 58, 495-508. Reinsch, C . H . (1967). Smoothing by spline functions. Numerische Mathematik, 10, 177-183. Shepp, L . A . (1966). Radon-Nykodym derivatives of Gaussian measures. Annals of Math-ematical Statistics, 37, 321-354. Shiller, R . J . (1984). Smoothness, priors and nonlinear regression. Journal of the American Statistical Association, 79, 609-615. Shilov, G . E . (1971). Linear Algebra. Prentice Hall , Englewood Cliffs, N J . Silverman, B . W . (1985). Some aspects of the spline smoothing approach to non-parametric regression curve fitting (with discussion). Journal of the Royal Statistical Society, Ser. B, 47, 1-52. Tong, H . (1990). Non-Linear Time Series: a Dynamical System Approach. Oxford sta-tistical science series, 6, Clarendon Press, Oxford. Wahba, G . (1978). Improper priors, spline smoothing and the problem of guarding against model errors in regression. Journal of the Royal Statistical Society, Ser. B, 40, 364-372. Wahba, G . (1983). Bayesian "confidence intervals" for the cross-validated smoothing spline. Journal of the Royal Statistical Society, Ser. B, 45, 133-150. 122 Wahba, G . (1990). Spline Models for Observational Data. S I A M , Philadelphia. Watson, M . W . and Engle, R, (1983). Alternative algorithms for the estimation of a dy-namic factor, M I M I C and varying coefficient regression models. Journal of Econometrics, 23, 385-400. Wecker, W . E . and Ansley, C . F . (1983). The signal extraction approach to nonlinear regression and spline smoothing. Journal of the American Statistical Association, 78, 81-89. Weinert, H . ; Byrd , R., and Sidhu, G . (1980). A stochastic framework for recursive com-putation of spline functions: part II, smoothing splines. Journal of Optimization Theory and Applications, 30, 2, 255-268. West, M . and Harrison, J . (1989). Bayesian Forecasting and Dynamic Models. Springer Series in Statistics, Springer Verlag, New York. W i l d , C . (1994). Embracing the "wider view" of statistics. The American Statistician, 48, 2, 163-171. Wold, H . O . (1938). A Study in the Analysis of Stationary Time Series, Almqvist & Wiksell, Uppsala, Sweden. 123 Appendix A Stochastic differential equations For a good review of this topic, see Arnold, 1974. Let (fi, T, P) be a probability space and W be an m-dimensional Wiener process (i.e. W = (W1,..., Wm)' where W* is Wiener process, i = 1 , . . . , n). Consider the filtration {!FS — o~(Wu, to < u < s)}, where cr(Q) represents the sigma algebra generated by the random variable Q. Define the space M;'m[t0,t] to be space of all d x m matrix-valued functions G, defined on [to,t] x 0, such that G(s, •) is ^-measurable Vs £ [to,t] and G(-,u) is with probability one in l ^ ^ o , t]-A . l Definition of the stochastic integral One wants to define [*GdW= [tG(s)dWs= f G{s,uj)dWs{u). Jto JtQ Jto First, the stochastic integral is defined for step functions G; i.e. there exists a partition to < ti < • • • < tn = t such that G(s) = (?(i;_i), Vs £ [ £ , _ ! , £ , ) , i = 1 , . . . , n . For such step functions, the stochastic integral of G with respect to Wt is the Rd random variable f GdW = YjG{ti_1){Wu-Wu_1). Jt° i=i In general, for any G function in M2'm[to,t] the stochastic integral of G with respect 124 to the m-dimensional Wiener process is defined as f GdW = P l im f GndW where {Gn} is a sequence of step functions in M ^ f ^ o , t] that approximates G in the sense of P l im / \G(s) - Gn(s)\2 ds = 0. p Here we denote X — P l im n—t-oo Xn if Xn • X as n > oo. A.2 Stochastic differentials The stochastic process defined by Xt(u>) = Xt0(u>) + f G{s,u)dWs{u), Jto Jto is said to possess the stochastic differential f(t) dt + G(t) dWt and this is denoted by writing dXt = f(t)dt + G(t)dW(t) = fdt + GdW. T h e o r e m A . 1 Let G denote a function in M2,m[to,T] and suppose that Xt= [tG{s)dWs, t0<t<T. Jt0 Then the following holds: • Xt is J-t-measurable. • If f E{\G{s)\2)ds < oo forallt<T, Jt0 then (XtjJ-'t), fort £ [to,?1], is an Rd-valued martingale; that is, for to < s < t < T, E(Xt\J-s) — Xt. Furthermore, fort,s € [io,?1] we have E{Xt) = 0, E{XtX's) = / E(G{u)G(u)') du; Jt0 in particular E(\Xt\2)= f E\G(u)\2 du. Jta • Xt has continuous sample path functions with probability 1. 125 A.3 Stochastic differential equations, existence and uniqueness of solutions Let Xt be an Rd stochastic process. A n equation of the form dXt = f(t, Xt) dt + G(t, Xt) dWu Xt0=c, t0<t<T <oo, (A.l) is called a stochastic differential equation. The random variable c is called the initial value at the instant to- Note that ( A . l ) is a symbolic way of writing Xt = c+ ff(s,Xs)ds + fG(s,Xs)dWs, t0<t<T<oo. (A.2) Jto Jto Remark A . l A n nth-order differential equation of the form dXt = f(t, Yt,..., y/""1)) + G(t, Yt,..., Y}n-'% with initial values Y^ = ct-, for i = 0 , 1 , . . . , n — 1, where Yt is i? d-valued and is an m-dimensional white noise, can be converted to the type ( A . l ) in the usual manner dXt = d Yt' Yl' \ I dt + [f(t,Yt,...,Y!n-V) 0 0 \ {G(t,Yt,...,Yt^) ) into a first-order stochastic differential equation for the i? d n -valued process Xt with initial value XtQ = ( co , . . . , c „ _ i ) ' . A stochastic process Xt is a solution of ( A . l ) if i) Xt is .^-measurable for t 6 [^o??1]; ii) the functions f(t,u) = f(t,Xt(u>)) and G(t,u) = G(t,Xt(u)) are such that with probability one / |/(s,w)| ds < oo and / \G(s,co)\ ds < oo\ Jto Jto iii) Equation (A.2) holds for every t G [to,T] with probability 1. 126 T h e o r e m A . 2 (Existence and uniqueness of the solution of a stochastic differential equa-tion) Consider the stochastic differential equation (A.1) where Wt is an Rm-valued Wiener process and c is a random variable independent of Wt — Wt0, for t > to- Suppose that the Rd-valued function f(t,x) and the d x ra matrix-valued function G(t,x) are defined and measurable on [to, T] x Rd and have the following properties: there exists a constant K > 0 such that: a) (Lipschitz condition) for all t € [t0, T], x G Rd, y £ Rd \f(t,x) - f(t,y)\ + \G{t,x) - G(t,y)\ < K\x - y\, where \G\2 = trGG'; b) (Restriction on growth) For all t € [to,T] and x € Rd, \f{t,x)\2 + \G{t,x)\2 <K2{\ + \x\2). Then, equation (A.1) has on [to,T] a unique Rd-valued solution Xt, continuous with prob-ability 1, that satisfies the initial condition Xto = c. If Yt is another solution of (A.1) then P ( sup \Xt-Yt\ > 0) = 0. \t0<t<T J A.4 Linear stochastic differential equations The stochastic equation (A.1) is said to be linear if f(t,x) = A(t)x + a(t), where A(t) is d x d matrix-valued and a(t) is Pcd-valued, and if G(t, x) = (B\[t)x + b\(i),..., Bm(t)x + bm(t)), where Bk(t) is d x d matrix valued and bk(t) is i? r f-valued. If a(t) — b\(i) — • • • = bm[t) = 0 then the linear stochastic differential equation is said to be homogeneous. If Bi(t) = • • • = Bm(t) = 0 the stochastic differential equation is said to be linear in the narrow sense. T h e o r e m A . 3 The linear (in the narrow sense) stochastic differential equation dXt = (A(t)Xt + a(t))dt + B(t)dWt, Xt0=c, (A.3) has on [to, T] the solution Xt = (c + f SOO"1^) ds + f §{s)-lB(s) dWs) ; (A.4) V JtQ Jt0 J where $ ( t ) is the fundamental matrix of the deterministic equation dXt/dt = A[t)Xt. 127 C o r o l l a r y 1 If the matrix A(t) = A in (A.3), then $(<) = e ^ ' - ' o ) = £ ~ 0 An(t-t0)n/n\ and X t = e ^ - ' o ) c + f* eA{t-s) [a,s) d s + B r s ) dWs}. (A>5) T h e o r e m A . 4 For the solution Xt of the linear stochastic differential equation dXt = (A(t)Xt + a(t)) dt + B{t) dWt, Xt0 = c, it holds that, under the assumption E(\c\2) < oo, • mt — E(Xt) = $(t) \E(C) + flo <&{s)~xa{s) ds . Therefore, mt is the solution of the deterministic linear differential equation dmt/dt = A{t)mt + a(t), mt0 = E(c). • The covariance function is K(s,t) = E[(X, - E(Xs))(Xt - E(Xt))'} rmin(s,t) = $ ( s ) Var(c)+ / ^(u^B^Biuymu))-1 du Jtn In particular, the variance matrix of the components of Xt, K{t) = K{t,t), is the unique symmetric nonnegative definite solution of the matrix equation dK(t)/dt = A{t)K(t) + K{t)A{t)' + B(t)B(t)', (A.6) with the initial value K(t0) = Var(c). T h e o r e m A . 5 The solution of the equation dXt = (A{t)Xt + a{t)) dt + B{t) dWt, Xt0 = c, is a stationary Gaussian process if A{t) = A, a(t) = 0, B{t) = B, the eigenvalues of A have negative real parts, and c is j\f(0, K)- distributed, where K is the solution /•oo K = / e BB'E 1 dt (A.7) Jo of the equation AK + KA' = —BB'. Then, for the process Xt, E(Xt) = 0 and cm(x,,x,) = E(x,,xl) = {e*yi\ ;|<|£ Under the above conditions, the process Xt is stationary in the wide sense with the above first and second moments even when c is not normally distributed but E(c) = 0 and E(cc') = K. 128 Appendix B Kronecker products Let A, B, C and D be matrices. The operator vec over the space oi n x m matrices into Rnm is defined as / a n • • • d m l \ vec \ a n l — • • • ? ®nli 021, . . . , Cl2ni • • • > ^ m l ) . . . , 0,mn)'; that is, the vec of a matrix is the operator that stacks the columns of the matrix into vector. The Kronecker product of A and B, denoted A® B, is defined as the matrix anP> • • • aimP> A®B a\nB B where A has entries a t j . Some useful properties of the vec operator and the Kronecker product are: • vec(o;A + B) = avec(A) + vec(B), for any a £ R. • vec(ABC) = (C ® A)vec(B). • (A®B)' = (A'®B'). • (A + B)®(C + D) = A®C + A®D + B®C + B®D. 129 • If A and B are invertible then A ® B is invertible and (A ® B)'1 = A - 1 ® B • (A® B)(C ®D) = (AC) ® (BD) and hence (A ® Bf - A2 ® B2. 130 Appendix C Jordan canonical form of a matrix Let j[k^ be a complex matrix denned by I A 1 0 A V 0 0 1 o \ 0 A 1 0 A J ( C l ) A matrix of the form ( C l ) is called a Jordan block. T h e o r e m C l Let A be an n x n complex matrix. Then, there are unique numbers A i , A 2 , . . . , \r, where A,- is an eigenvalue of A, i = l , . . . , r , and unique positive integers m i , m 2 , . . . , mr such that A is similar to the matrix J = diag(jl?\jJrK..,jM), (C.2) that is, there exist a nonsingular matrix P such that P~lAP = J. The form (C.2) is called the Jordan canonical form of the matrix A. It is unique up to the ordering of the Jordan blocks. Note that the same eigenvalue can appear in different blocks. 131 Appendix D The matrix exponential For good references on this topic, see for example Golub and Van Loan (1983) and the references therein. Given a matrix A , the matrix exponential of A is formally defined as oo eAt = J2(At)k/k\. The matrix eAt is always defined and invertible with ( e A t ) _ 1 = e~At. If A is invertible then f etA dt = A _ 1 e t j 4 . Another property of the matrix exponential is that if the matrices A and B commute then eA+B = eAeB = eBeA. T h e o r e m D . l Let A = P J P - 1 be a n x n matrix with Jordan canonical form J as in (C.2). Then, eAt = Pdiag(eJi\...,eJrt)p-\ where J{ = J>™'\ i = 1,... ,r. Example: Let A = d i a g ( d a , . . . ,dn). Then, etA = d i a g ( e W l , . . . ,etdn) Example: Let E be the n x n matrix given by 0 1 E 1 0 132 Then one can easily show that E , k <r, where Ek has l's on the kth superdiagonal and is zero elsewhere. Ek = 0 if k > r. Hence, 1 t t2/2 . . . tn-ll{n-l)\ t2/2 t 1 Example: Let A be a Jordan block of the form ( C . l ) . Then, A = XI + E, where E is as in the previous example and 1 t t2/2 ... i " - V ( n - l ) ! 0 •• '•• '•• : ; ••. ••. ••. t2/2 : '•• '•• t 0 0 1 „At Xt e — e Example: Let A 0 co - u 0 Then eAt = cos(u;tf) sin(a;t) — sin(u;i) cos(ujt) • Some software that is useful for computing matrix exponentials is Maple and Octave. Maple computes matrix exponentials in a symbolic way. Octave is free software and for more information, the reader should visit the Internet sites by http://www.che.wise, edu/octave/readme.html and http://csugrad. cs. vt. edu/manuals/octave/octave-10.html. The author of Octave is John W . Eaton, Department of Chemical Engineering, University of Wisconsin. 133 Appendix E Data sets Data from Draper and Smith (1981, p.476) on fraction of available chlorine of a product. Length of time Average since produced Available Available (Weeks) Chlorine Chlorine 8 0.49,0.49 0.490 10 0.48,0.47,0.48,0.47 0.475 12 0.46,0.46,0.45,0.43 0.450 14 0.45,0.43,0.43 0.437 16 0.44,0.43,0.43 0.433 18 0.46,0.45 0.455 20 0.42,0.42,0.43 0.423 22 0.41,0.41,0.40 0.407 24 0.42,0.40,0.40 0.407 26 0.41,0.40,0.41 0.407 28 0.41,0.40 0.405 30 0.40,0.40,0.38 0.393 32 0.41,0.40 0.405 34 0.40 0.400 36 0.41,0.38 0.395 38 0.40,0.40 0.400 40 0.39 0.390 42 0.39 0.390 134 Data from Andrews and Herzberg (1985) on the age-adjusted incidences of mel the years 1936-1972. Year Incidence Year Incidence per 100,000 per 100,000 1936 0.9 1955 2.9 1937 0.8 1956 2.5 1938 0.8 1957 2.6 1939 1.3 1958 3.2 1940 1.4 1959 3.8 1941 1.2 1960 4.2 1942 1.7 1961 3.9 1943 1.8 1962 3.7 1944 1.6 1963 3.3 1945 1.5 1964 3.7 1946 1.5 1965 3.9 1947 2.0 1966 4.1 1948 2.5 1967 3.8 1949 2.7 1968 4.7 1950 2.9 1969 4.4 1951 2.5 1970 4.8 1952 3.1 1971 4.8 1953 2.4 1972 4.8 1954 2.2 135
Thesis/Dissertation
eng
Statistics
Vancouver : University of British Columbia Library
For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use.
Graduate
The model based approach to smoothing
Text