KALMAN FILTER AND ITS APPLICATION TO FLOW FORECASTING by PATRICIA NGAN B. A. Sc., Uni v e r s i t y of B r i t i s h Columbia, 1981 M. S., C a l i f o r n i a I n s t i t u t e of Technology, 19&2 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY i n FACULTY OF GRADUATE STUDIES ( C i v i l Engineering) We accept t h i s thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA A p r i l 1985 P a t r i c i a Ngan, 1985 In presenting t h i s thesis i n p a r t i a l f u l f i l m e n t of the requirements for an advanced degree at the University of B r i t i s h Columbia, I agree that the Library s h a l l make i t f r e e l y available for reference and study. I further agree that permission for extensive copying of t h i s thesis for scholarly purposes may be granted by the head of my department or by h i s or her representatives. I t i s understood that copying or publication of t h i s thesis for f i n a n c i a l gain s h a l l not be allowed without my written permission. Department of 6lQvt- &Nl6t^G:eftlM| The University of B r i t i s h Columbia 1956 Main Mall Vancouver, Canada V6T 1Y3 5 - 6 f ^ / R - n ABSTRACT The Kalman F i l t e r has been applied to many f i e l d s of hydrology, p a r t i c u l a r l y in the area of flood forecasting. This recursive estimation technique i s based on a state-space approach which combines model description of a process with data information, and accounts for uncertainties in a hydrologic system. This thesis deals with applications of the Kalman F i l t e r to ARMAX models in the context of streamflow prediction. Implementation of the Kalman F i l t e r requires s p e c i f i c a t i o n of the noise covariances (Q, R) and i n i t i a l conditions of the state vector (x 0, P 0 ) . D i f f i c u l t i e s arise in streamflow applications because these quantities are often not known. Forecasting performance of the Kalman F i l t e r i s examined using synthetic flow data, generated with chosen values for the i n i t i a l state vector and the noise covariances. An ARMAX model is cast into state-space form with the c o e f f i c i e n t s as the state vector. S e n s i t i v i t y of the flow forecasts to s p e c i f i c a t i o n of x 0, P 0, Q, R, (which may be d i f f e r e n t from the generation values) i s examined. The f i l t e r ' s forecasting performance i s mainly affected by the combined s p e c i f i c a t i o n of Q and R. When both noise covariances are unknown, they should be s p e c i f i e d r e l a t i v e l y large in order to achieve a reasonable forecasting performance. S p e c i f i f y i n g Q too small and R too large should be avoided as i t results in poor flow forecasts. i i The f i l t e r ' s performance i s also examined using actual flow data from a large r i v e r , whose behavior changes slowly with time. Three simple ARMAX models are used for t h i s investigation. Although there are di f f e r e n t ways of writing the ARMAX model in state-space form, i t is found that the best forecasting scheme i s to model the ARMAX c o e f f i c i e n t s as the state vector. Under t h i s formulation, the Kalman F i l t e r i s used to give recursive estimates of the c o e f f i c i e n t s . Hence flow predictions can be revised at each time step with the la t e s t state estimate. This formulation also has the feature that i n i t i a l values of the ARMAX co e f f i c i e n t s need not be known accurately. The noise variances of each of the three models are estimated by the method of maximum l i k e l i h o o d , whereby the li k e l i h o o d function i s evaluated in terms of the innovations. Analyses of flow data for the stations considered in this thesis, indicate that the variance of the measurement error i s proportional to the square of the flow. In practice, flow predictions several time steps in advance are often required. For autoregressive processes, t h i s involves unknown elements in the system matrix H of the Kalman model. The Kalman algorithm underestimates the variance of the forecast error i f H and x are both unknown. For the AR(1) model, a general expression for the mean square error of the forecast i s developed. It i s shown that the formula reduces to the Kalman equation for the case where the system matrix i s known. The importance of t h i s formula i s realized in forecasting situations where management decisions depend on the r e l i a b i l i t y of flow predict i o n s , reflected by their mean square errors. iv Table of Contents Chapter Page ABSTRACT i i LIST OF TABLES x LIST OF FIGURES x i i ACKNOWLEDGEMENTS x i i i 1. INTRODUCTION 1 1.1 Problem d e f i n i t i o n 1 1.2 General Objectives 2 1.3 General Thesis Outline 6 2. LITERATURE REVIEW 8 2.1 Explanation of the Kalman F i l t e r 8 2.2 General Applications 12 2.3 Hydrologic Applications 13 2.4 Problems in hydrologic applications 14 2.4.1 Unknown dynamics 14 2.4.2 Estimation of State system matrix 15 2.4.3 Unknown Observation system matrix 16 2.4.4 Estimation of the noise covariances 16 2.5 Detailed thesis objectives and outline 19 3. SENSITIVITY ANALYSIS 22 3.1 Experimental Plan 23 3.2 Experimental Procedure 25 3.2.1 Generation of Flow Data 25 3.2.2 Kalman F i l t e r Specifications 27 3.2.3 F a c t o r i a l Experiment 28 3.2.4 Analysis of Variance 32 v 3.3 -Results 34 3.3.1 ANOVA Results 34 3.3.2 Two-dimensional Graphs 36 3.3.3 Three-dimensional plot 42 3.4 Discussion of Results 44 3.4.1 ANOVA Results 44 3.4.2 Interpretation of Graphs 45 3.5 Summary 45 4. HYDROLOGIC SYSTEMS 48 4.1 Time-invariant Systems 49 4.2 Three ARMAX Models 51 4.3 Scope of Applications 53 5. ESTIMATION OF MEASUREMENT NOISE VARIANCE 54 5.1 The Estimation Method 54 5.2 Estimation Procedure 56 5.3 Transformation of flow data to achieve Uniform Variance 57 5.3.1 Relation between the Noise Variance and the Flow 58 5.3.2 Results 59 5.3.3 Conclusions 60 5.3.4 Box-Cox Transforms 60 5.4 Maximum Likelihood Estimation of the Noise Var iance • 62 5.4.1 Theory 62 5.4.2 Results 64 5.5 Summary 65 6. APPLICATION: AR(1) 67 v i 6.1 Description of the Forecasting schemes 68 6.1.1 Properties of Scheme 1 68 6.1.2 Properties of Scheme 2 69 6.1.3 Properties of Scheme 3 69 6.2 General Discussion 72 6.3 Experimental Procedure 73 6.4 Results 75 6.4.1 Estimate of AR c o e f f i c i e n t 75 6.4.2 Estimate of the noise variance 76 6.4.3 Performance Indicators for Scheme 1 76 6.4.4 Performance Indicators for Scheme 2 78 6.4.5 Performance Indicators for Scheme 3 80 6.5 Discussion of Results. 83 6.5.1 Forecasting Performance of Scheme 1 83 6.5.2 Forecasting Performance of Scheme 2 84 6.5.3 Forecasting Performance of Scheme 3 84 6.6 Summary 87 7. APPLICATION: TRANSFER FUNCTION MODEL 88 7.1 Description of Forecasting Schemes 90 7.1.1 Properties of Scheme 1 90 7.1.2 Properties of Scheme 2 90 7.2 General Discussion 91 7.3 Experimental Procedure 92 7.4 Results 93 7.4.1 Estimate of the model c o e f f i c i e n t s 93 7.4.2 Estimate of the noise variance 94 7.4.3 Performance Indicators for Scheme 1 95 v i i 7.4.4 Performance Indicators for Scheme 2 97 7.5 Discussion of Results 99 7.6 Summary 100 APPLICATION: ARMAX MODEL 101 8.1 Properties of the Forecasting Scheme 102 8.2 General Discussion 103 8.3 Experimental Procedure 104 8.4 Results 105 8.4.1 Estimate of the noise variance 105 8.4.2 Performance indicators for the ARMAX model 106 8.4.3 Comparison of 1-step Performance Indicators 108 8.4.4 Comparison of 2-step Performance Indicators 109 8.5 Discussion of Results 111' 8.5.1 Performance of the ARMAX model 111 8.5.2 Comparison of 1-step performance 111 8.5.3 Comparison of 2-step performance 111 8 . 6 Summary 113 VARIANCE OF THE FORECAST ERROR 114 9.1 Background 115 9.2 Problem d e f i n i t i o n 115 9.3 Covariances of Products of normal random variables 117 9.4 Variance of the forecast error for the scalar case 118 9.5 I l l u s t r a t i o n of the variance formula for the AR(1) 120 9.6 Conclusions 122 vi i i 10. SUMMARY AND CONCLUSIONS 125 10.1 Summary of Thesis Contributions 126 10.2 S e n s i t i v i t y Analysis 129 10.3 Maximum Likelihood Estimation of the noise variance 130 10.4 The Autoregressive Model 130 10.5 Transfer Function Model 132 10.6 Combined ARMAX Model 133 10.7 Variance of the forecast error 134 10.8 State-space Formulation 135 10.9 Future Directions 135 BIBLIOGRAPHY 137 ix LIST OF TABLES Table No. T i t l e Page 3.1 Values of x 0*, Q*, R* used in flow 26 generation 3.2 Schematic Coding Scheme 27 3.3 Number of Levels for the Factors 28 3.4 Levels of the Noise Covariances 29 3.5 Computer Output - ANOVA Results for CODE 7 35 5.1 Values of the slope for years 1970-1979 60 5.2 Maximum Likelihood Estimates of the Noise 64 Var iance 5.3 Averaged Maximum Likelihood Estimate of the 66 Noise Variance 6.1 Properties of the forecasting schemes 71 6.2 Least Squares Estimates of the 75 autoregressive c o e f f i c i e n t 6.3 Maximum Likelihood Estimates of the noise 76 variance 6.4 Range of values for input parameters 77 6.5a Values of PI, for Scheme 1 77 6.5b Values of PI 2 for Scheme 1 77 6.5c Values of PI 3 for Scheme 1 78 6.6 Range of values for a 78 6.7a Values of PI, for Scheme 2 79 6.7b Values of PI 2 for Scheme 2 79 6.7c Values of PI 3 for Scheme 3 80 6.8 Range of values for input parameters 81 x 6.9a Values of PI, for Scheme 3 81 6.9b Values of P I 2 for Scheme 3 82 6.9c Values of PI 3 for Scheme 3 82 7.1 Least Squares Estimates of the model 94 co e f f i c i e n t 7.2 Maximum Likelihood Estimates of the noise 95 variance 7.3 Performance Indicators for Scheme 1 95 7.4 Performance Indicators for the AR(1) model 96 7.5 S e n s i t i v i t y of the Performance Indicators 96 to the model c o e f f i c i e n t 7.6 Performance Indicators for Scheme 2 97 7.7 Range of Values for input parameters 98 7.8 Performance Indicators for Scheme 2 98 8.1 Maximum Likelihood Estimates of the noise 105 var iance 8.2a Variation of PI, to noise s p e c i f i c a t i o n 106 8.2b Variation of PI 2 to noise s p e c i f i c a t i o n 106 8.2c Variation of PI 3 to noise s p e c i f i c a t i o n 107 8.3a Variation of PI , to noise s p e c i f i c a t i o n 107 8.3b Variation of PI 2 to noise s p e c i f i c a t i o n 107 8.3c Variation of P I 3 to noise s p e c i f i c a t i o n 108 8.4 Values of PI for ARMAX and AR(1) models 109 8.5 Values of PI for ARMAX and Transfer 110 Function models xi LIST OF FIGURES F i g . No. T i t l e Page 1.1 Geographic location of the study area. 4 1.2 Outflow hydrograph at Hope, B.C. 5 2.1 A Schematic representation of the Kalman 11 F i l t e r procedure at each time step. 2.2 S e n s i t i v i t y of the state estimates to 14 mis-specifications of noise standard deviations. (After Mehra, 1968) 3.1 Schematic diagram of Experimental 33 Procedure. 3.2a RRMS with respect to Q (CODES 1-4) 38 3.2b RRMS with respect to R (CODES 1-4) 38 3.3a RRMS with respect to Q (CODES 5-8) 39 3.3b RRMS with respect to R (CODES 5-8) 39 3.4a RRMS with respect to Q (CODES 9-10) 40 3.4b' RRMS with respect to R (CODES 9-10) 40 3.4a RRMS with respect to Q (CODES 11-12) 41 3.4b RRMS with respect to R (CODES 11-12) 41 3.5 Response Surface - RRMS with respect to QR 43 speci f i c a t i o n s 4.1 A 'data window' approach to s l i g h t l y 50 time-invariant hydrologic systems. x i i ACKNOWLEDGEMENTS I wish to thank Dr. S.O. (Denis) Russell (Department of C i v i l Engineering), and Dr. Piet de Jong (Faculty of Commerce) for their valuable guidance and advice throughout t h i s research. The knowledge I have gained i s deeply appreciated, and the experience has been most rewarding. Sincere thanks go to Dr. M.C. Quick (Department of C i v i l Engineering) and Dr. W. Caselton (Department of C i v i l Enginnering) for their comments on streamflow modelling during t h i s work. Appreciation i s given to Mr. R.G. Boals (Environment Canada, Water Resources Branch, Ottawa), and Mr. Oliver Nagy (Water Survey of Canada, Vancouver) for their prompt attention to the flow data r e q u i s i t i o n . Thanks are also extended to Dr. Glen Cooper (Computing Centre at UBC) and Mr. Gordon Finlay (Department of C i v i l Engineering, Graphics Lab) for their assistance in the computing aspects of t h i s research. I am grateful to my friends for their informative comments throughtout t h i s thesis. F i n a l l y I wish to thank my mother for her support and encouragement, without which th i s work would not have been possible. xi i i 1. INTRODUCTION Flow forecasting is an important aspect in the operation and control of water resource systems. S t a t i s t i c a l models are used extensively in the stochastic modelling of riverflows. Adaptive forecasting occurs when predictions can be updated at each time step in response to incoming observations. The Kalman F i l t e r i s an example of an adaptive forecasting scheme. This estimation method i s based on a linear state-space model. Because time series or regression models are amenable to state-space formulation, they are preferred over conceptual models such as SWMM, HEC I, and HEC II in adaptive hydrologic forecasting. Applications of the Kalman F i l t e r have been very successful in communications and aerospace engineering, f i e l d s in which the system dynamics and the governing physical equations are well known. However, such is not the case in streamflow applications. The performance of the Kalman F i l t e r may be greatly affected when these equations and dynamics are unknown. 1.1 Problem d e f i n i t i o n Although the Kalman F i l t e r i s explained in Chapter 2, i t i s b r i e f l y defined here in order to state the general problem addressed in t h i s thesis. The linear state-space model, sometimes known as the Gauss-Markov model consists of two equations: State eqn: xfc = • * t 21 -^1 + - t wfc^(0,Q) 1.1 1 2 Observation egn: y_t = Hfc x_t + Y_t v / n/(0,R) 1.2 The objective i s to recursively estimate the state vector, xfc, based on the current observations, y_t. The Kalman algorithm provides a method whereby state estimates are continually updated as new observations are received at each time step. The state vector i s corrupted by white noise, w.tA>(0,Q) with covariance matrix Q. Simil a r l y , the vector of observations i s corrupted by white noise v^/(0,R) with covariance matrix R. The algorithm assumes that H, Q, and R are known at time t. A f i l t e r e d estimate x i s given by the Kalman F i l t e r at every time step. In addition, the mean square error for the state estimate is given. A problem which often arises in practice i s that of unknown dynamics; either system matrices H or noise covariances Q , R. As O'Connell and Clarke (1981) point out, this d i f f i c u l t y i s often overlooked in applications even though the state estimation procedure depends c r i t i c a l l y on the proper s p e c i f i c a t i o n of dynamics. 1.2 G e n e r a l O b j e c t i v e s This research work deals with the problem of unknown dynamics. ARMAX models for describing streamflow phenomenon are considered in t h i s thesis. These are time series models with exogeneous inputs. 3 The p r a c t i c a l application used in t h i s thesis i s motivated by flood forecasting in B.C. The chosen problem i s streamflow prediction for the Fraser River at Hope. This station i s located near the downstream portion of the Fraser River, (see F i g . 1.1) As a r e s u l t of the large drainage area (217,000 km2), flow magnitudes are in the order of thousands m3/s. Peak flow usually occurs in June, with a t y p i c a l value of 5000 m3/s. F i g . 1.2 gives an outflow hydrograph for 1983. S t a t i o n n e a r Spences B r i d g e S t a t i o n at Texas Creek /Thompson \ \ R i v e r / W. Vancouver' .TOLL CAMtOlA ISIA*0' "A* ""VANCOUVER ^ O Richmond n VAIOU 'O 1 i n = 4 5 km = 28 mi F i g . 1.1 G e o g r a p h i c l o c a t i o n o f the s t u d y a r e a . OUTFLOW HYDROGRAPH 8000 I Time F i g . 1. 2 O u t f l o w h y d r o g r a p h a t Hope, B. 6 In t h i s t h e s is, ARMAX models are cast into the state-space format with the c o e f f i c i e n t s as the state vector. In hydrologic applications, future flow i s the quantity of interest. Forecasting performance of the f i l t e r i s measured in terms of the observation forecast error (y-y). The general objectives of the thesis are: 1. to examine the s e n s i t i v i t y of the Kalman F i l t e r with respect to mis-specifications of noise covariances and i n i t i a l conditions of the state vector. These inputs are required by the Kalman algorithm, but are often unknown in p r a c t i c e . 2 . to investigate the maximum l i k e l i h o o d method for estimating the noise covariances in hydrologic models of p r a c t i c a l interest. 3. to compare the forecasting performance of ARMAX models depending on whether the Kalman F i l t e r i s used or not. 1.3 General Thesis Outline The next chapter reviews the l i t e r a t u r e on the Kalman F i l t e r . Included in this review are applications and the problems which arise in practice. Chapter 3 investigates the s e n s i t i v i t y of the f i l t e r with respect to input quantities, which are assumed known by the algorithm. Chapter 4 considers three special cases of the ARMAX models which are most often used in streamflow predictions. In subsequent chapters, the research deals with these three models. Chapter 5 describes the maximum l i k e l i h o o d method for 7 estimating the noise variance. Chapters 6, 7, and 8 examine the forecasting performance of the three models, depending on whether or not the Kalman F i l t e r i s used. Often in forecasting, the standard deviation of the prediction i s required; but there are many cases where the expression for the variance given by the Kalman F i l t e r i s not applicable. Chapter 9 discusses t h i s in d e t a i l and a t h e o r e t i c a l expression i s presented for the variance of the forecast error. The importance of t h i s i s i l l u s t r a t e d in an example. F i n a l l y , chapter 10 summarizes the results of t h i s thesis and gives recommendations for future research. 2. LITERATURE REVIEW There are fiv e sections to t h i s chapter. F i r s t , the Kalman F i l t e r i s described in d e t a i l . General applications are given in section 2 and hydrologic examples in section 3. P r a c t i c a l problems in the implementation of the f i l t e r are discussed in section 4. F i n a l l y , section 5 outlines in d e t a i l the s p e c i f i c topics addressed in subsequent chapters. 2.1 Explanation of the Kalman F i l t e r The following notation i s used in t h i s t h e s i s . Vectors are denoted by underlined, small l e t t e r s ; and matrices are denoted by c a p i t a l l e t t e r s . The linear state-space model, sometimes known as the Gauss-Markov model consists of two equations: State eqn: xfc = $ X f i + w w /vR(0,Q) 2.1 The elements of xfc are the state variables to be estimated, but they need not be physically measurable. These states are modelled as a Markov process. The manner in which they evolve through time i s given by 3>, the state t r a n s i t i o n matrix. The states are corrupted by white noise, w d i s t r i b u t e d with mean 0 and covariance matrix Q. $ and Q can be time-varying, but are assumed known at a l l times. Observation eqn: Y_fc = H xfc + y_t v^H{0 ,R) 2.2 8 9 Measurements taken at time t are related to the states of the system through H. They are also corrupted by white noise, y_t d i s t r i b u t e d with mean 0 and covariance matrix R. Sim i l a r l y , H and R may be time-varying, but are assumed known at a l l times. Moreover, wfc and v f c are assumed to be s e r i a l l y and mutually uncorrelated. For the Kalman model consisting of eqn. 2.1 and 2.2, the objective i s to estimate the state vector at current time, xfc given a l l past measurements. Obtaining an estimate for the state at the current time, i s known as f i l t e r i n g . In pa r t i c u l a r , i f a l l past information can be summarized in a prior estimate for the state vector, t h i s type of estimation i s known as recursive f i l t e r i n g . The Kalman F i l t e r i s a recursive procedure which yiel d s estimates for xfc at every time step. This f i l t e r i n g technique can be given a Bayesian interpretation. Posterior estimates of the state variables are obtained by updating their prior estimates through the measurements received at time t. However, i t i s not a Bayesian technique per se, as often quoted in the l i t e r a t u r e , despite t h i s decision theoretic approach. In fact, x t (the 'posterior' or updated estimate) i s derived from a least squares c r i t e r i o n by minimizing the residual of the state estimation vector. The Kalman F i l t e r i s a solution which yie l d s a l i n e a r , minimum variance estimator for the state vector at every time step. In addition, the error covariance matrix P , associated with x., i s given. 10 The Kalman F i l t e r algorithm can be divided into two pa r t s: 1 . Prior to observing y_ -t/t-1 = *t-1 -t-1 2 , 3 P t / t - 1 = *t-1 Pt-1 *'t-1 + Q t 2 ' 4 2. After observing y_t, i t = h/t-i + K t [ * t - £ t ] 2 - 5 P t = [I - Kfc H t] 2.6 where Kfc = P ^ ^ H't i ^ t / t ^ » \ + * t ]' 1 2.7 The 1-step forecast for the observation i s : £t+1 = H A + l / t 2 ' 8 The Kalman F i l t e r accounts for uncertainties by providing an algorithm which sequentially combines model (state equation) and data (measurement equation) information to y i e l d updated estimates of the state vector. These estimates can be projected forward to obtain future predictions of the observations. It i s computationally appealing due to i t s recursive nature, i . e . a l l previous information i s contained in the prior estimate of the state. A schematic representation of th i s procedure i s shown in F i g . 2.1. 11 S t a t e e q u a t i o n at time t-1 1 ~ Q O b s e r v a t i o n s J • at time t R i F i g . 2.1 A Schematic r e p r e s e n t a t i o n of the Kalman F i l t e r procedure at each time s t e p . 12 2.2 General Applications Application of the Kalman F i l t e r i s widespread, and the l i t e r a t u r e spans over many d i s c i p l i n e s . It i s used p a r t i c u l a r l y in forecasting situations. The c l a s s i c example i s from aerospace engineering. The problem i s to continually estimate the position, v e l o c i t y and acceleration of a target in Cartesian co-ordinates, upon receiving noise corrupted measurements from a radar in polar co-ordinates. In finance, Kalman F i l t e r i n g has been applied to estimate the regression c o e f f i c i e n t s in a model for stock earnings per share (Mehra, 1979). This technique has also been applied in water quality physical models involving BOD-DO equations. Examples of t h i s are found in the work by Young and Whitehead (1977), Constable and McBean (1979). However, no n - l i n e a r i t i e s in these formulations lead to the use of the extended Kalman F i l t e r which i s not the subject of the thesis. The Kalman F i l t e r has also been i l l u s t r a t e d in many aspects of the hydrologic process. Bras (1978) has successfully applied this technique to sampling network design. An example of estimating groundwater basin c h a r a c t e r i s t i c s i s given by McLaughlin (1978). However, the dimensionality of groundwater models i s large, as many parameters need to be estimated. Computational e f f i c i e n c y and the problem of joint state and parameter estimation are two major d i f f i c u l t i e s (Wilson et a l . , 1978). This f i l t e r i n g technique has also been used in the operation of water resource systems. An example i s the 13 maintaining of water l e v e l s for navigation purposes. Duong et al.(l978) applied the Kalman F i l t e r to obtain real time estimates of system parameters required for c o n t r o l l i n g a lock and dam gate. 2.3 Hydrologic Applications The most abundant applications of thi s estimation technique are in streamflow forecasting. In many situations, the underlying process can be described by a s t a t i s t i c a l model such as a regression or time series model. In these cases, the Kalman F i l t e r i s a useful supplement to hydrologic forecasting. The r a i n f a l l - r u n o f f process i s often modelled by an autoregressive scheme plus some input information. There are two ways of formulating t h i s type of hydrologic model into the state-space notation. Rodriguez-Iturbe (1978) has formulated the problem with the c o e f f i c i e n t s of an ARMA model as the variables of interest. This i s also shown in an example by Wood et a l . , (1978). The instantaneous unit hydrograph, IUH represented by a convolution i n t e g r a l , i s a c l a s s i c example of such a forecasting scheme. For examples, see Rodriguez-Iturbe (1978), Szollosi-Nagy (1976). A l t e r n a t i v e l y , the discharges can be taken as the states, but the unknown ARMAX parameters then need to be estimated a p r i o r i . An example of t h i s i s given in Szollosi-Nagy et al. , ( l 9 7 7 ) . Conceptual response models involving the continuity and storage discharge equations have also been 14 used with the Kalman F i l t e r . However, these models lead to non-linear estimation. Wood (1978) investigated flood routing models via the Kalman F i l t e r for forecasting water l e v e l s . 2 .4 Problems i n h y d r o l o g i c a p p l i c a t i o n s 2 .4.1 Unknown dynamics The problem of unknown noise covariances i s often encountered in hydrologic modelling as a result of the complexity of the runoff process. Mehra (1979) pointed out the extreme s e n s i t i v i t y of the Kalman F i l t e r to underspecification of R, or the misspecification of Q. 100. in = 2 I oc . — E Ui IU 0.1 Forward .Smoothing I 0.5 10 0, asumad/o. actual Effect of errors i n measurement noise variance on steady-state RMS error i n state e s t i -mation. or-= measurement noise standard deviation, (from Mehra 1968). J _ attumad/O. actuat q q Effect of errors i n process noise variance on steady-state RMS error in state estimation. o°n = process noise standard deviation, (from Mehra 1968). F i g . 2.2 S e n s i t i v i t y of the state estimates to mis - spec i f i ca t ions of noise standard dev iat ions . 15 This was shown in q u a l i t a t i v e terms and the idea conforms with one's i n t u i t i o n . However, the f i l t e r ' s performance under simultaneous misspecification of Q and R needs to be examined. In practice, these noise covariances are spec i f i e d too small or too large with respect to their true values. The manner in which these s p e c i f i c a t i o n s a f f e c t the f i l t e r ' s performance i s important knowledge for the forecaster. Moreover, the i n i t i a l conditions of the state are often unknown. These questions are addressed in the next chapter. Because flow i s the variable of interest, the f i l t e r ' s performance indicators are based on the observation forecast error. 2 . 4 . 2 Estimation of State system matrix Various methods have been proposed in the past for estimating these unknown matrices, H, Q, and R. $ i s the state system matrix. When elements of # are unknown, thi s leads to non-linear f i l t e r i n g . One way of resolving t h i s i s to l i n e a r i z e the problem with respect to the lat e s t state estimate by the Extended Kalman F i l t e r . However, Mehra has indicated that the estimates thus obtained are sensitive to i n i t i a l conditions. Moreover, the EKF has not been proven to y i e l d consistent estimates (Young, 1977). In an example on the Ombrone Basin in Italy, Szollosi-Nagy (1977) used the instrumental variable (IV) approach by Young (1977) to resolve the no n - l i n e a r i t i e s of the ARMA c o e f f i c i e n t s in An extension to t h i s approach, developed by Todini et 16 al.,(l978) known as MISP technique was also applied to the Ombrone Basin. This made use of two Kalman F i l t e r s : one which y i e l d s the state estimates (given the parameter values) and the other which yie l d s parameter estimates (given the state values). A constrained linear estimator proposed by Natale and Todini (1976) was i l l u s t r a t e d in an application by Wood (1978). 2.4.3 Unknown Observation system matrix H is the system matrix in the observation equation of the Kalman model. When elements of H are unknown, thi s necessarily leads to a di f f e r e n t expression for the variance of the forecast error. Authors have overlooked t h i s (Harrison & Stevens, 1976) when they suggested using the same Kalman F i l t e r formula to obtain the variance. Whether or not thi s i s a serious error in practice depends on the applicati o n . Feldstein (1971) has studied t h i s problem in the context of econometrics. However, because one of his assumptions i s often not s a t i s f i e d in hydrologic applications, this problem needs further investigation. 2.4.4 Estimation of the noise covariances Implementation of the Kalman F i l t e r allows the forecaster to exercise his judgement regarding the accuracy of the underlying model vs. that of the observations. This i s achieved by specifying the noise covariances, Q and R at each time step. The Kalman gain matrix K, in e f f e c t , 17 controls how the state estimates are updated. Hence the performance of the f i l t e r depends c r i t i c a l l y on these values. These are often unknown in water resource applications, and many authors have addressed t h i s problem. In real-time applications where h i s t o r i c a l data are scarce, adaptive algorithms are often used (Mehra 1973, IEEE). Jazwinski & B a i l i e (1969), Sage & Husa (1969) propose covariance matching techniques, where the t h e o r e t i c a l covariance of the innovations and the sample covariance are matched by adjusting Q and R appropriately. This h e u r i s t i c approach i s computationally a t t r a c t i v e though i t i s not guaranteed to converge. The special case of known Q but unknown R i s reported to have been handled more successfully (Mehra, 1972). Correlation methods only apply to time-invariant systems under stationary conditions which impose a l i m i t a t i o n for some water resource systems. The Bayesian method i s used to compute the posterior p r o b a b i l i t y of the 6 = {Q, R} given the observed data, i . e . p(0|y f c). Calculation of t h i s requires the l i k e l i h o o d of 6, £(fl|yt) = p(y t| y t_.j r . . .0). If there are p elements in 6,and each has N choices, then there are N*3 possible combinations of the set 6 corresponding to Kalman F i l t e r s . Posterior p r o b a b i l i t y i s calculated for each one of these Kalman F i l t e r s . As i t can be readily seen, th i s i s not a t t r a c t i v e computationally unless N and p are both small. This procedure discriminates the best model out of a given group 18 of models. However, i f the o r i g i n a l set of 0's do not encompass the true value, then the "best model" i s s t i l l not optimal in the sense that there i s s t i l l a better 6 that could be used. This method was used by Moore & Jones (1978), and Valdes et a l . , (1978). The tremendous complexity of the system renders i t computationally unattractive and cumbersome. The previously described methods are somewhat ad hoc in the i r applications. The Maximum l i k e l i h o o d method for estimating Q and R i s based on a mathematical p r i n c i p l e and can be universally applied. The main advantage of th i s method i s that i t i s asymptotically unbiased, consistent and e f f i c i e n t . The idea i s to maximize the l i k e l i h o o d of observing a p a r t i c u l a r combination of Q and R given the observation y over a range of values for the noise covariances. The major drawback i s that the l i k e l i h o o d expression involves computation with a non-diagonal matrix. Because the previously described methods are not guaranteed to work, t h i s method i s investigated in chapter 5. It should be noted that although the problem of estimating noise variances has been addressed t h e o r e t i c a l l y in the l i t e r a t u r e , the results have had li m i t e d use in pr a c t i c e . The main reasons are: 1. Many of the proposed solutions were developed for some par t i c u l a r application. Hence they tend to be application dependent. 2. The methods which are e a s i l y implemented computationally 19 are usually h e u r i s t i c in nature. They have not been proven to give r e l i a b l e estimates. Moreover, d i f f e r e n t authors, who have used them, reported contradictory results (Mehra, 1972; O'Connell & Clarke, 1981; Sage & Husa, 1969; Wood et^ a l . , 1978). 3. The more complex methods are not e a s i l y implemented, though in general, they y i e l d more consistent estimates (Mehra, 1972,1980). in addition, the complex mathematics required do not lend themselves to p r a c t i c a l implementations. In fact, they act more l i k e a black-box to the forecaster. There i s a d e f i n i t e need to bridge the gap between the theoretical research and the implementation f e a s i b i l i t y required in practice. 2 .5 Detailed thesis objectives and outline As discussed in section 2.4, the noise covariances and the i n i t i a l s p e c i f i c a t i o n s for the state vector are often unknown. These correspond to Q, R, x 0, Po of the Kalman state-space model. Two questions are addressed in chapter 3. 1. Which of the above quantities (Q, R, x 0, P 0) i f misspecified, have a p r a c t i c a l e f f e c t on the f i l t e r ' s performance. 2. How do these quantities a f f e c t the forecasting performance. The second half of the thesis i s based on ARMAX models. S t a t i s t i c a l models are chosen for the applications because 20 they are amenable to state-space formulation of the Kalman F i l t e r . Conceptual models such as SWMM and HEC I are not considered in this thesis. These are physically based models and are more d i f f i c u l t to implement. These are more worthwhile to consider for long-term forecasting for a river basin. On the other hand, ARMAX models require r e l a t i v e l y less data. This thesis i s focused on the Kalman F i l t e r rather than model i d e n t i f i c a t i o n of hydrologic processes. Thus, s t a t i s t i c a l models are used. Moreover, the simpler representation of streamflow phenomenon allows a more direc t investigation of the p r a c t i c a l problems associated with the Kalman F i l t e r . Chapter 4 describes three ARMAX models to be used for the remainder of the thesis. The method of maximum li k e l i h o o d i s used to estimate the noise variance for these models in chapter 5. In p a r t i c u l a r , a s i m p l i f i e d approach to evaluating the log- l i k e l i h o o d expression i s examined. This has been proposed by Schweppe and i l l u s t r a t e d by Ledolter ett a l . (1983). The next three chapters investigate the forecasting performance of the three models depending on whether or not the Kalman F i l t e r i s used. Chapter 6 considers an AR(1) model for predicting flows 1 day ahead at Hope. Chapter 7 examines a transfer function model for predicting the 2-day advance flow. Chapter 8 studies a combined model which gives 1 and 2 day forecasts. The 1 and 2 step forecasts are compared to those obtained from the models of chapters 6 and 7 respectively. F i n a l l y , a general 21 expression for the mean square error of the forecast when H is unknown, i s developed in Chapter 9. This s i t u a t i o n occurs, for example, when the AR(1) model i s used to predict the flow than more two time steps ahead. In t h i s case, the variance of the forecast error i s not given by the Kalman equation, HPH' + R. 3. SENSITIVITY ANALYSIS The Kalman F i l t e r has been used in many hydrologic applications. It i s a popular estimation technique due to i t s recursive nature and i t s a b i l i t y to handle uncertainties. However, i t s applications have not always been successful. This i s often because the assumption of known input quantities, i s not s a t i s f i e d . These quantities are the i n i t i a l conditions of the state vector x 0 , P 0, and the noise covariance matrices Q, R. The Kalman algorithm assumes that x 0t P o i Q» R are c o r r e c t l y s p e c i f i e d . But when they are unknown in practice, their mis-specification may result in unreliable estimates. The objective of t h i s study i s to examine the s e n s i t i v i t y of the f i l t e r to mis-specification of the quantities x 0, P 0, Q, R. The s e n s i t i v i t y study i s analyzed as a s t a t i s t i c a l f a c t o r i a l experiment. Each of the input quantities {x 0, P 0, Q, R} i s treated as a factor. The forecasting performance of the f i l t e r i s measured by an indicator based on the observation forecast error. This indicator i s the response variable in the f a c t o r i a l exper iment. A scalar ARMAX model i s used for streamflow modelling in t h i s study. In t h i s application, an autoregressive model of order 1, with temperature as an exogeneous variable i s chosen. The model i s cast into state-space format with the model c o e f f i c i e n t s as the state vector. Hence Q represents the noise covariance matrix for the c o e f f i c i e n t s , and R i s 22 23 the variance of the measurement error. Results of t h i s s e n s i t i v i t y study show that i n i t i a l s p e c i f i c a t i o n s of the ARMAX c o e f f i c i e n t s , x 0 and P 0, do not aff e c t the forecasting performance of the f i l t e r to any degree of p r a c t i c a l concern. The performance i s materially affected by the combined s p e c i f i c a t i o n of the noise covariances, Q and R. In addition, conclusions are made as to how the unknown noise covariances should be specified in order to achieve good forecasting performance. Specifying both Q and R r e l a t i v e l y large results in a f i l t e r performance comparable to the case of known Q and R. Specifying Q small and R large y i e l d s the worst forecasting performance. If the measurement noise variance i s known, then the forecasting performance i s in d i f f e r e n t to under or over-specification of Q. However, i f Q is known, i t i s found that better f i l t e r performance i s obtained i f the measurement noise variance i s sp e c i f i e d r e l a t i v e l y large. 3.1 Experimental Plan The approach to t h i s investigation i s divided into four parts. 1. Flow generation Flow sequences are generated by a chosen ARMAX model with known input quantities: x 0*, Q*, R*. 2. Kalman F i l t e r S p e c i f i c a t i o n The flow sequences are treated as observations and a Kalman F i l t e r i s applied to each flow sequence. The 24 state-space format used corresponds to allowing the ARMAX c o e f f i c i e n t s be the state vector in the Kalman model. Hence, estimates of the c o e f f i c i e n t s are given at each time step and flow forecasts are made with the latest estimate given by the f i l t e r . The noise covariance matrices used in the generation of streamflow data are Q* and R*. The Kalman algorithm uses x 0, P 0, Q, and R where these may be d i f f e r e n t from the generation values. Q and R can be s p e c i f i e d to be less than, equal to, or greater than their true values, Q* and R*. Hence, the s e n s i t i v i t y of the f i l t e r ' s performance to these spe c i f i c a t i o n s i s examined. The performance of the f i l t e r i s measured by an indicator based on the observation forecast error. A Kalman F i l t e r i s applied to each combination of the input quantities. One of these f i l t e r s contain the correct s p e c i f i c a t i o n s of the inputs. F a c t o r i a l Experiment The investigation i s studied as a s t a t i s t i c a l f a c t o r i a l experiment with four factors; namely {x 0, P 0, Q, R}. The levels of these factors correspond to the d i f f e r e n t s p e c i f i c a t i o n s . Analysis of Variance S e n s i t i v i t y of the f i l t e r ' s performance i s investigated through analysis of variance (ANOVA) on the performance indicators. The analyses indicate which factors or combinations thereof, have a s i g n i f i c a n t effect on the 25 performance indicator. 3.2 Experimental Procedure 3.2.1 Generation of Flow Data Application of the Kalman F i l t e r in t h i s thesis is motivated by real-time flow forecasting of the Fraser River. A stochastic model for describing the streamflow phenomenon is chosen: q t = a t *t-1 + b t T e m p t + v t 3.1 The autoregressive term i s c h a r a c t e r i s t i c of stations with large drainage areas, as the flow i s dependent on upstream storage e f f e c t s . The temperature term represents the snowmelt influence on the runoff. The c o e f f i c i e n t s of the model, a f c and b f c, are chosen to follow a random walk. The ARMAX model of eqn. 3.1 i s recast into state-space framework as follows: State eqn: V " a t - r + "w, t" b t _ bt-1_ _ w 2 t . 3.2 Obs. eqn: [ q t _ 1 Tempt] 3.3 The noise terms are chosen with the following properties: 26 1• wfc ~ (0,Q*) The covariance matrix for the c o e f f i c i e n t s (states) is constant with respect to time. It is also a diagonal matrix. Thus, independent errors for the model c o e f f i c i e n t s are assumed. 2. v t<v (0,R*) S i m i l a r l y , a constant variance for the observation error is used. Batches of flow sequences are generated for each testing condition denoted as "CODE". Each sequence has T=100 time steps. Each CODE corresponds to a pa r t i c u l a r combination of x_o*, Q* i R* • The following table gives the values used in the flow generation. Table 3.1 Values of x 0*, Q*, R* used in flow generation x 0 * = .85 7.00 Arbitrary s t a r t i n g values are used as t h i s does not affe c t the streamflow generation. Q i * = 000036 0 0 .0016 0 * = y I I .0001 0 0 .01 0 * = y I I I .0004 0 0 .0225 R:* = 100 R n * = 400 = 625 R I V* = 900 The chosen values of Q* are representative of p r a c t i c a l s i t u a t i o n s ; a standard deviation of .5% and 2% from one time 27 step to the next. S i m i l a r l y , the range of R* r e f l e c t s a measurement error with a standard deviation from 3% to 12% of the flow. In practice, the standard deviation of the measurement error i s not more than 10% of the flow. The number of combinations of the above input quantities i s twelve. Hence, the number of CODES i s also twelve. Table 3.2 gives a schematic layout of the coding scheme. Table 3.2 Schematic coding scheme V R * II R * III R * IV CODE 1 CODE 2 CODE 3 CODE 4 0 * y I I CODE 5 CODE 6 CODE 7 CODE 8 0 * U I I I CODE 9 CODE 10 CODE 11 CODE 12 3.2.2 Kalman F i l t e r Specifications Values of the quantities used in data generation are denoted with a '*'; otherwise the l e t t e r s represent specifications of these quantities for the Kalman algorithm. Therefore, Q* i s used in generation; while Q i s a sp e c i f i c a t i o n value. Different s p e c i f i c a t i o n s for {x 0i Por Qf R l result in d i f f e r e n t Kalman F i l t e r s . For each 28 code, 36 Kalman F i l t e r s are evaluated, each with a sample size of 5 streamflow sequences. The number of streamflow sequences required under each CODE is 180. For each of the 180 series, the performance of the Kalman F i l t e r i s measured in terms of the flow forecast error. The indicator i s rel a t i v e root mean square (RRMS) error defined as: This i s preferred over the common RMS error as RRMS does not depend on the units of measurement. 3.2.3 F a c t o r i a l Experiment The investigation i s set up as a f a c t o r i a l experiment with a completely randomized design. It is a fully-crossed four factor experiment. The factors are the spe c i f i c a t i o n s of the input quantities: {x0r po> Q> R) • The number of levels for each factor i s given below: Table 3.3 Number of Levels for the Factors RRMS 3.4 Q R 2 2 3 3 l 29 The l e v e l s correspond to d i f f e r e n t s p e c i f i c a t i o n s of the input quantities. The number of ways of combining these levels i s : k = II / . = 36. Hence, 36 d i f f e r e n t Kalman i F i l t e r s are used under each CODE. The levels of these factors are described below. Table 3.4 Levels of the Noise Covariances le v e l x 0 P 0 Q, R 1 g< ?od es ~o. 9 r _6.50_ >timate conf idei guess ".001 0" _0 . 065_ it too small .25Q*, .25R* 2 b< ad esti "0.425-1 0.50 .mate non-con guess \75 0~ Eident correct spec i f icat ion Q*, R* 3 too large 4Q*, 4R* The good estimate of x 0 corresponds to an i n i t i a l s p e c i f i c a t i o n error of 7%, while the bad estimate corresponds to 50%. These l i m i t s are representative of the range of accuracies for estimates of x 0 in practi c e . Level 1 of the noise s p e c i f i c a t i o n corresponds to under-estimating the standard deviation by a factor of 2. Si m i l a r l y , l e v e l 3 30 corresponds to over-estimating by a factor of 2. Level 2 represents a s p e c i f i c a t i o n which i s equal to the true covariance matrix. In practice, i t i s expected that estimates of Q and R, are within the bounds provided by lev e l s 1 and 3. The response variable of the f a c t o r i a l experiment i s the performance indicator, RRMS. Its values are analyzed by ANOVA to order to determine which of the four factors or combinations thereof, have an important e f f e c t on the performance of the Kalman F i l t e r . Each of the 36 Kalman F i l t e r s i s subjected to n=5 series of observations, the 5 series being d i f f e r e n t in each case. This conforms to a completely randomized (CR) design. For t h i s non-repeated measures design, the measurement errors are uncorrelated. Thus, an unrepresentative series does not re p l i c a t e i t s c h a r a c t e r i s t i c from one Kalman F i l t e r to another. A l t e r n a t i v e l y , the same 5 series of observations could be used for each Kalman F i l t e r in order to minimize the v a r i a t i o n between groups. This repeated measures design results in correlated measurement errors among the various Kalman F i l t e r s . Hence, a CR design i s chosen. The s t a t i s t i c a l model for a f a c t o r i a l experiment with 4 factors i s : Y. i j klm = M + A. + B.+C, + D i l k + 2-way interactions .. + 3-way interactions ... + 4-way interaction + e m 3.5 31 For t h i s s e n s i t i v i t y study, the following variable substitution i s made: RRMS = M + X 0 + P 0 + Q + R x_oPo + • • • + x 0P 0Q + ... + x 0P 0QR + e The notation i s described as follows: 1. RRMS i s the response variable of the experiment. 2. ju i s the o v e r a l l mean of the performance indicator, RRMS. 3. x 0, P 0, Q, R are the factors and they are c a l l e d the main effects in eqn. 3.6 4. The 2-way interactions include a l l possible combinations of the factors in groups of two. These are x 0P 0, x 0Q, x 0R, P0Q, P0R, QR. 5. S i m i l a r l y , the 3-way interactions are a l l combinations of the factors in groups of three. These are x 0P 0Q, x 0P 0R, P 0QR, x 0QR. 6. F i n a l l y , there i s the 4-way interaction and the error, 3.6 32 3.2.4 Analysis of Variance The experiment is analyzed by ANOVA with a l l terms (except M) on the right hand side of eqn. 3.6 as the 'sources' in the ANOVA table. One analysis i s done for each CODE; hence, twelve ANOVA's are made in t o t a l . The computer package 'ANOVAR' i s used and a sample output of the ANOVA table i s given in the Results (Section 3.3.1). The output gives the p-value for each source under 'F PROB'. These ANOVA results are used as a q u a l i t a t i v e guide for pointing out those sources which are important in p r a c t i c a l situations. The p-values are compared among the sources in order to indicate those which have a s i g n i f i c a n t effect on the response, RRMS. The smaller the p-value, the more l i k e l y i t i s that the corresponding source has a s i g n i f i c a n t effect on the f i l t e r ' s performance. A flow chart of the procedure i s given in F i g . 3.1. 33 G e n e r a t e f l o w s choose x R'" o g e n e r a t e N=kn=180 f l o w s e q u e n c e s o f 100 time s t e p s l o n g i Kalman F i l t e r . K F 1 , 1 K F 1 , 2 •• • K F 1 , 5 S p e c i f i c a t i o n s K F 3 6 , 1 K F 3 6 , 2 " • K F 3 6 , 5 F a c t o r i a l E x p e r i m e n t measure the p e r f o r m a n c e of the Kalman F i l t e r , RRMS ANOVA d e t e r m i n e s i g n i f i c a n t s o u r c e s f o r RRMS. Repeat f o r each CODE F i g . 3 . 1 S c h e m a t i c d i a g r a m o f E x p e r i m e n t a l P r o c e d u r e . 34 3.3 Results 3 . 3 . 1 ANOVA Results The ANOVA results indicate that under a l l CODES, the si g n i f i c a n t sources are the s p e c i f i c a t i o n of the main effect Q, and the interaction QR. The significance of QR means that i t i s the combined s p e c i f i c a t i o n of the noise covariances which a f f e c t s the f i l t e r ' s performance. A sample output of the ANOVA table for CODE 7 i s shown on the following page. The p-values (represented by F PROB in the table) are compared. Those corresponding to the sources Q and QR are an order of magnitude lower than the rest. This phenomenon is noted under a l l CODES. The manner in which the QR interaction a f f e c t s RRMS i s obtained from the following two and three-dimensional graphs. The mean RRMS values are plotted with respect to the leve l s of the source, QR. 35 T a b l e 3.5 Computer O u t p u t ANOVA R e s u l t s f o r CODE 7 F a c t o r s : A = * 0 B = P 0 C = Q D = R p - v a l u e s ; g i v e n by F PROB ANALYSIS OP VARIANCE/COVARIANCE FOR VARIABLE RRMS SOURCE A B AB C AC BC ABC D AD BO ABD CD ACD BCD ABCD ERROR TOTAL SUM OF D.F. SOUARES 1 1.439339E-05 1 3.582272E-05 1 2.683472E-05 2 3.062712E-03 2 5.325221E-04 2 \ 1.691381E-04 2 8.595134E-04 2 9.816524E-04 2 8.214938E-04 2 5.839004E-04 2 5.987364E-04 4 4.261532E-03 4 2.407318E-03 4 1.477050E-03 4 1.577517E-03 144 4.805392E-02 179 6.54S426E-02 MEAN SOUARE 1.439339E-05 3.582271E-05 2.683472E-05 1.531356E-03 2.662609E-04 S'.456905E-05 4.2975S7E-04 4.9092G2E-04 4.1074S8E-04 2.919501E-04 2.993S81E-04 1.O65383E-03 6.018295E-04 3.692624E-04 3.943790E-04 3.337075E-O4 F VALUE F PROB 0. 0431 0. 8177 0. 1073 0. 7398 0. 0804 0. 7689 4 . 5689 0. 01 17 0. 7979 0. 4559 0. 2534 0. 7791 1 . 2878 0. 2785 1 . 471 1 0. 2316 1 . 2309 0. 2949 0. 8749 0. 4219 0. 8971 0. 4126 3. 1926 0. 0151 1 . 8035 0. 1300 1 . 10S5 0. 3560 1 . 1818 0. 3211 Note t h a t t h e p - v a l u e s f o r s o u r c e s C=Q and CD=QR a r e .0117 and .0151 r e s p e c t i v e l y . 36 3.3.2 Two-dimensional Graphs The ordinate in the following graphs i s the performance indicator, RRMS. This indicator i s such that the smaller the value, the better i s the performance. The abscissa denotes the three levels of s p e c i f i c a t i o n for the noise covariance. For instance, Q(1) corresponds to specifying Q too small. Graphs in series (a) are plotted with respect to Q for each l e v e l of R ( i . e . holding R constant). Those in series (b) are with respect to R for each l e v e l of Q. The features displayed by the graphs are summarized below: 1. From series (b) of Figs. 3.2-3.4, i t can be seen that the smallest RRMS value occurs at Q(2)R(2). This represents correct s p e c i f i c a t i o n s of the noise covariances. 2. The graphs in series (a) show that i f R i s given c o r r e c t l y , then specifying Q too small Q(1), or too large Q(3), results in almost i d e n t i c a l RRMS values. In other words, the f i l t e r ' s performance i s ind i f f e r e n t to mis-specification of Q i f R i s known. 3. A d i f f e r e n t behavior i s observed by examining the graphs in series (b) of Figs. 3.2-3.4 for le v e l Q(2). It i s noted that i f Q i s given c o r r e c t l y , then under-specification of R results in a larger RRMS than over-specification. This phenomenon i s p a r t i c u l a r l y noticeable under Q* I I I shown in F i g . 3.4, for systems with large disturbances in the states. 4. For each CODE, nine RRMS values are plotted. In each 37 case, two of these values are r e l a t i v e l y close to the minimum RRMS. These values correspond to the sp e c i f i c a t i o n s , Q(2)R(3) and Q(3)R(3). It i s noted that both these s p e c i f i c a t i o n s involve R(3). This feature i s more pronounced as Q* and R* increase. In general, specifying R too small should be avoided. High values of the performance indicator for a l l s p e c i f i c a t i o n s of Q under l e v e l R(1) are obtained. A large increase in the RRMS value is noted for Q(1)R(3). For example, see F i g . 3.3 (b). 38 F i g . 3.2 ( a ) RRMS w i t h r e s p e c t t o Q F i g . 3.2 (b) RRMS w i t h r e s p e c t - t o R 39 CODE 5 CODE 5 0.03*. 0.03-t Legend 0.08 CODE 8 Legend o ot » Q2_ & 03 F i g . 3.3 ( a ) RRMS w i t h r e s p e c t to Q F i g . 3.3 (b) RRMS w i t h r e s p e c t t o R 40 F i g . 3.4 ( a ) RRMS w i t h r e s p e c t t o Q F i g . 3.4 (b) RRMS w i t h r e s p e c t t o R F i g . 3.4 ( a ) c o n t ' d RRMS w i t h r e s p e c t t o Q F i g . 3.4 (b) c o n t ' d RRMS w i t h r e s p e c t t o R 42 3.3.3 Three-dimensional plot A three-dimensional plot of the interpolated response surface for the RRMS i s given. This i s plotted with respect to the nine s p e c i f i c a t i o n levels of both Q and R. Because the behavior of the f i l t e r to QR interaction i s the same, only one response surface i s drawn. This i s a q u a l i t a t i v e plot as the axes are scaled appropriately to give a reasonable image. The surface minimum occurs in the region corresponding to Q(2)R(2). The response surface r i s e s in a l l directions away from the minimum. It i s noted that the topography i s not too d i f f e r e n t from the minimum, in a zone defined by Q(2)R(3), and Q(3)R(3). The surface also shows a sharp increase in RRMS for the s p e c i f i c a t i o n , Q(1)R(3). Q l e v e l s RRMS R l e v e l s A3 Q ( 1 ) R ( 3 ) F i g . 3.5 Response S u r f a c e RRMS w i t h r e s p e c t t o QR s p e c i f i c a t i o n s 44 3.4 Discussion of Results 3.4.1 ANOVA Results The f a c t o r i a l experiment considers four factors: x 0, P 0, Q, R. It i s found that: 1. I n i t i a l s p e c i f i c a t i o n s of the ARMAX c o e f f i c i e n t s x 0, P 0 do not have a s i g n i f i c a n t effect on the f i l t e r ' s forecasting performance. In practice, these c o e f f i c i e n t s are often unknown and t r a d i t i o n a l l y , much work has gone into their estimation p r i o r to forecasting. This study shows that even an error of 50% in t h e i r i n i t i a l s p e c i f i c a t i o n , does not s i g n i f i c a n t l y degrade the performance. 2. With the c o e f f i c i e n t s as the state variables in the Kalman model, i t i s found that the f i l t e r i s sensitive to the combined s p e c i f i c a t i o n of the noise covariances, Q and R. The 1-step flow forecast i s based on the l a t e s t estimate of the state vector, which is" updated by using the Kalman gain matrix, Kfc. Gelb (1979) shows that the Kalman gain can be interpreted as a r a t i o of the uncertainties in the state estimate to the measurement noise. Because the uncertainty in the state estimate i s affected by the noise covariance Q, Kfc can be considered as a r a t i o of Q to R. Hence, i t i s the combined s p e c i f i c a t i o n of the noise covariances which affects the f i l t e r ' s performance. 45 3.4.2 Interpretat ion of Graphs Because of the interaction between Q and R, i t i s not s u f f i c i e n t to examine the f i l t e r ' s s e n s i t i v i t y to either Q or R i n d i v i d u a l l y . It should be examined with respect to both qua n t i t i e s . 1. The optimum le v e l of performance (represented by the minimum RRMS value) i s obtained when Q and R are spe c i f i e d c o r r e c t l y . This is true under a l l testing conditions. 2. If the variance of the measurement error R i s known, then the f i l t e r ' s performance i s i n d i f f e r e n t to under or overspecification of Q. Its performance level can only be improved i f Q i s sp e c i f i e d c o r r e c t l y . 3. If Q i s known to begin with, then specifying R too small should be avoided. 4. For systems with large disturbances in the ARMAX c o e f f i c i e n t s , under-specification of Q should be avoided. This i s p a r t i c u l a r l y true when R i s sp e c i f i e d large. The worst performance occurs for Q(1)R(3). 5. A general guideline i s to specify both covariance matrices large because r e l a t i v e l y good forecasting performance can be obtained under Q(3)R(3). 3.5 Summary The Kalman F i l t e r i s applied to an ARMAX model for streamflow forecasting. A common problem i s that i n i t i a l conditions of the model c o e f f i c i e n t s and the noise 46 c o v a r i a n c e s are unknown. The s e n s i t i v i t y of the f i l t e r ' s f o r e c a s t i n g performance to the s p e c i f i c a t i o n of x 0 , P 0 , Q, R i s examined. The ARMAX model i s r e c a s t i n t o s t a t e - s p a c e format wi th the c o e f f i c i e n t s as the s t a t e s i n the Kalman model . For t h i s s e n s i t i v i t y s t u d y , the system m a t r i c e s $ and H of the Kalman F i l t e r are known. The f i l t e r ' s performance i s measured in terms of the o b s e r v a t i o n f o r e c a s t e r r o r , d e f i n e d by RRMS. The - c o n d i t i o n s t e s t e d c o r r e s p o n d to a s t a n d a r d d e v i a t i o n of .5 to 2% change i n the model c o e f f i c i e n t s , and a s t a n d a r d d e v i a t i o n of 3 to 12% i n the measurement n o i s e . I t i s found tha t i n i t i a l s p e c i f i c a t i o n s of the ARMAX c o e f f i c i e n t s do not a f f e c t the f i l t e r ' s per formance . T h i s f e a t u r e i s u s e f u l f or the p r a c t i c i n g e n g i n e e r ; e s p e c i a l l y when t h e r e i s i n s u f f i c i e n t data f o r o b t a i n i n g r e l i a b l e e s t i m a t e s of the c o e f f i c i e n t s p r i o r to f o r e c a s t i n g . The s tudy a l s o shows that the f i l t e r ' s performance i s s e n s i t i v e to the combined s p e c i f i c a t i o n of the n o i s e c o v a r i a n c e s . I f Q and R are unknown, i t i s best to s p e c i f y both c o v a r i a n c e s l a r g e r than t h e i r expected v a l u e s . T h i s s p e c i f i c a t i o n r e s u l t s i n a f o r e c a s t i n g performance comparable to the case of known Q and R. I f o n l y one of the n o i s e c o v a r i a n c e s can be e s t i m a t e d , i t i s b e t t e r to e s t imate R a c c u r a t e l y . T h i s i s due to the f i l t e r ' s i n s e n s i t i v i t y to under or over s p e c i f i c a t i o n of Q, i f R i s g i v e n c o r r e c t l y . However, i f Q i s g i v e n c o r r e c t l y , the f i l t e r performs worse i f R i s s p e c i f i e d too s m a l l . 47 In general, i t i s best to specify a large variance for the measurement error. This results in good forecasting performance as long as Q i s not specified much smaller than i t s true value. 4. HYDROLOGIC SYSTEMS The Kalman F i l t e r has been applied in the previous s e n s i t i v i t y study to a general ARMAX model, used to generate synthetic streamflow data. The state-space formulation corresponds to l e t t i n g the ARMAX c o e f f i c i e n t s be the state variables in the Kalman model. The study shows how the forecasting performance of the f i l t e r i s affected by the combined noise covariance s p e c i f i c a t i o n . It also indicates that i n i t i a l s p e c i f i c a t i o n of the ARMAX c o e f f i c i e n t s do not af f e c t the flow forecasting performance. This chapter describes three p a r t i c u l a r ARMAX models, chosen for Kalman F i l t e r applications in Chapters 6 to 8. Actual d a i l y flow data i s used in these studies for forecasting flows 1 and 2 days in advance. These models represent hydrologic systems with certain basin c h a r a c t e r i s t i c s , described in thi s chapter. The p r a c t i c a l i t y of the Kalman F i l t e r cannot be explored for a l l types of hydrologic systems represented by the general ARMAX model, as i t i s data dependent to some extent. Hence, s p e c i a l i z a t i o n to certain types of systems i s necessary. The p r a c t i c a l i t y of the Kalman F i l t e r and i t s performance under d i f f e r e n t state-space formulations of a hydrologic model are investigated using real data. Studies in Chapters 6 to 8 use d a i l y flow data from the Fraser River at Hope in B.C. The Fraser River at Hope i s one of the largest rivers in Canada and the flow changes r e l a t i v e l y slowly through time. As discussed in section 4.1, 48 49 th i s allows the assumption of time-invariant systems. Therefore, ARMAX models with constant c o e f f i c i e n t s w i l l be considered for these studies. In practice then, only the noise variance associated with the measurement equation in the Kalman model need to be estimated prior to forecasting. This subject is addressed in Chapter 5. The three ARMAX models chosen are r e l a t i v e l y simple from a hydrologic point of view because they do not involve complicated physical equations. Physically based models have been used for some types of r i v e r basins (e.g. Quick and Pipes, 1976). These are usually deterministic models requiring extensive data inputs. In addition, the system of equations do not readily lend themselves to state-space formulation. Because the focus of subsequent studies i s on the Kalman F i l t e r and i t s potential in flow forecasting rather than on developing the best forecasting procedure, time-invariant ARMAX models are selected. These represent the simplest types of hydrologic systems which are of p r a c t i c a l interest in water resource management. 4.1 Time-invariant Systems Flow forecasting using ARMAX models with time-invariant c o e f f i c i e n t s are considered. Constant c o e f f i c i e n t models are appropriate for hydrologic systems whose behavior remains constant with time. These types of systems are usually characterized by large drainage areas, where the flows are influenced by storage e f f e c t s . For these systems, r a i n f a l l inputs have l i t t l e influence on the outflow. Rodriguez-Iturbe (1978) and Szollosi-Nagy (1976) concluded that the assumption of time-invariant c o e f f i c i e n t s can be extended to systems whose response c h a r a c t e r i s t i c s change slowly over time. A "window" type of approach i s suggested for handling the time varying nature, as shown in the diagram. TIME ( f r o m S z o l l o s i - N a g y , 1976) z ( t ) i s the o b s e r v a t i o n at time t . F i g . 4'. 1 A 'd a t a window' a p p r o a c h t o s l i g h t l y t i m e - v a r i a n t h y d r o l o g i c s y s t e m s . 51 Application of the Kalman F i l t e r in this thesis i s motivated by the following problem. Streamflow predictions for the Fraser River at Hope are required in the control of thi s water resource system. The maximum dail y flow of the year (peak runoff) which occurs in June, i s mainly due to snowmelt. Therefore, a reasonable forecasting horizon i s A p r i l 1 to September 30. Within t h i s time frame, i t i s expected that the basin c h a r a c t e r i s t i c s do not change abruptly. The station at Hope i s located near the downstream end of the Fraser River. This means that i t has a large drainage area and thus a slow basin response. Hence, time-invariant ARMAX models are appropriate for the investigations in thi s thesis. As a result, only one noise variance i s to be estimated, and thi s i s addressed in Chapter 5. 4.2 Three ARMAX Models Three stochastic models commonly used in streamflow modelling are considered. 1. AR(1) q t = aq t_ 1 + v f c 2. TRANSFER FUNCTION q t = b i q * t - 2 + b2<3**t-2 + V t 3. Combined, ARMAX q t = c 1 q t _ 1 + c 2 q * t _ 2 + c 3q** t_2 + v f c where q, q*, q** denote streamflows at three d i f f e r e n t 52 stations. These are special cases of the general ARMAX model. They can be thought of as a (p,q) ARMA model with (r) exogeneous inputs. General j u s t i f i c a t i o n s for using these types of models have been given in the previous section. Particular reasons for using each of these models are discussed below. 1. The autoregressive model of order 1 has been widely used in streamflow modelling. Quimpo (1973) has shown that storage e f f e c t s can be represented by autoregressive terms. Predictions beyond t h i s can be made i f the unknown flows at lag one are replaced by their estimate. However, th i s necessarily results in a larger variance for the flow forecast error. This topic is addressed in Chapter 9. Application of the Kalman F i l t e r to this model i s investigated in Chapter 6. 2. Upstream inputs can contribute to the r i s i n g portion of an outflow hydrograph as pointed out by Wood (1978). A transfer function model can be used where the co e f f i c i e n t s b, and b 2 represent flow contributions from upstream stations. The 2-day lag in the flows allow the prediction of the 2-day forecast with a l l the exogeneous variables known at time t. The application of the Kalman F i l t e r to thi s model i s examined in Chapter 7. In fact, flood routing procedures described by the convolution integral are often l i n e a r i z e d . The discrete time representation of the int e g r a l can be considered as a transfer function model. An example of th i s 53 formulation can be found in Wood(l978). 3. The t h i r d model combines the autoregressive process with upstream inputs. The Kalman F i l t e r i s applied to this model in Chapter 8. 4.3 Scope of A p p l i c a t i o n s The state-space formulation whereby the ARMAX c o e f f i c i e n t s form the elements of the state vector i s considered. This format i s used throughout the remaining chapters, except where otherwise noted. The problems associated with the s p e c i f i c a t i o n of the inputs {x 0, P 0, Q, R} have been reduced to the problem of estimating R only because: 1. Q=0 as time-invariant c o e f f i c i e n t s are considered. 2. I n i t i a l s p e c i f i c a t i o n of the ARMAX c o e f f i c i e n t s , x 0, P 0 do not affe c t the f i l t e r ' s forecasting performance. 5. ESTIMATION OF MEASUREMENT NOISE VARIANCE Results of Chapter 3 show that the f i l t e r ' s forecasting performance i s affected by the s p e c i f i c a t i o n of the noise covariance matrices. In fact, variance estimation in state-space models i s a subject of much current i n t e r e s t . A review of t h i s was given in Chapter 2, section 2.4.4. In this chapter, a method for estimating the noise variance i s investigated. The three models discussed in Chapter 4 are considered. It i s found that the variance of the measurement error cannot be assumed to be constant with time. Analysis indicates that i t varies in proportion to the square of the flow. Therefore, a Box-Cox transform is used for s t a b i l i z i n g the variance. The appropriate transformation i s to take logarithms of the o r i g i n a l measurements. The new series has a noise variance which is approximately constant. The method of maximum l i k e l i h o o d i s then used to estimate the noise variance for each model. An equivalent l o g - l i k l i h o o d function for the flow i s written in terms of the observation forecast error; and i s evaluated through the use of the Kalman F i l t e r . 5.1 The Estimation Method Because an accurate and r e l i a b l e estimate of the noise variance i s desired, h e u r i s t i c methods discussed in Chapter 2 are not considered. More sophisticated methods e x i s t , but are t r a d i t i o n a l l y avoided due to extensive computational 54 55 requirements. However, the a v a i l a b i l i t y of computer software, has made these methods more a t t r a c t i v e . An example is the Bayesian method of estimation. This w i l l not be used in t h i s thesis but i s discussed b r i e f l y as i t has attracted some attention in the recent l i t e r a t u r e . This method does not estimate the noise variance per se. Actually, i t i s a discrimination technique. Numerous Kalman models are considered at the s t a r t ; each one containing a possible value of the noise variance, a2^. It i s assumed that one of these Kalman F i l t e r s contains the correct noise variance a*2. Posterior p r o b a b i l i t y of a2 ^ being o*2 given the observations, i s calculated for each model. After s u f f i c i e n t number of observations, the model with a2- closest to a*2 ' 1 w i l l be assigned the highest posterior p r o b a b i l i t y . Although i n i t i a l applications by Valdes et a l . , (1978) show that this method does select the model with a2 ^ closest to the true variance, i t is not considered here for the following reasons: 1. Applications of th i s technique i s only in i t s i n i t i a l stage, hence i t s properties s t i l l need to be examined. 2. In a study by the previous authors, i t was found that 300 time steps are required before the best model can be distinguished. Since we are only interested in forecasting for part of the year, there is not enough data for the discrimination procedure. 3. The o r i g i n a l set of models i s assumed to encompass the true value of R. 56 This l a s t point presents a l i m i t a t i o n in practi c e . Often, the engineer is not s u f f i c i e n t l y familiar with the physical c h a r a c t e r i s t i c s of a basin to provide the bound values. Hence, the set of a 2^'s chosen may not include the true value. Instead, the method of maximum l i k e l i h o o d i s used in th i s chapter to estimate the noise variance. The maximum l i k e l i h o o d method gives estimates which are known to be consistent and asymptotically e f f i c i e n t . No prior knowledge of the system parameters are required here. The d i f f i c u l t y with t h i s method i s the evaluation of the li k e l i h o o d function. In time series applications, errors in successive streamflow measurements are l i k e l y to be correlated. Harvey (1981) outlines in his text, how the evaluation of the log l i k e l i h o o d function can be si m p l i f i e d through the use of the Kalman F i l t e r . An example of thi s estimation procedure i s i l l u s t r a t e d by Ledolter (1983) in a Management Science context. 5.2 Estimation Procedure An engineer often has to make decisions for a water resource system in real time. It i s desirable that the variance of the measurements errors be considered constant with respect to time. This can be approximately achieved by transforming the raw measurements. For each of the models which w i l l be used in subsequent studies, the noise variance i s estimated by the following procedure. 1. Determine the appropriate Box-Cox transformation for the 57 o r i g i n a l flow data such that the noise variance of the new series i s approximately constant. 2. For the transformed data, obtain the maximum l i k e l i h o o d estimate of the noise variance by using the Kalman F i l t e r . 5.3 Transformation of flow data to achieve Uniform Variance Simple applications of stochastic flow models assume that the measurement errors are constant for the raw flow data. For an example, see Wood (1978). However, a more r e a l i s t i c assumption i s that the variance of the flow observation error be a function of the flow i t s e l f . Patry and Marino (1983) examined various noise assumptions commonly used in hydrologic models. In an application by Natale & Todini (1976), a linear dependence of the noise term's standard deviation on the flow i s assumed. In a real-time forecasting s i t u a t i o n , flow prediction i s made simpler i f the measurement errors are time-invariant. Transforms for s t a b i l i z i n g the variance are avai l a b l e . A Box-Cox transform is used in thi s application. In general, th i s class of transforms i s appropriate when the variance of the observations can be expressed as a function of i t s mean. Before the appropriate Box-Cox transformation can be determined, the r e l a t i o n between the variance of the observations and i t s mean i s investigated. 58 5.3.1 Relation between the Noise Variance and the Flow The three ARMAX models of Chapter 4 are used in subsequent chapters for the application of the Kalman F i l t e r . The state-space formulation where the ARMAX co e f f i c i e n t s are the state variables are considered. The conditional variance of the flow, given the model c o e f f i c i e n t s , i s equal to the variance of the measurement error. This i s expressed as the following, where z represents the o r i g i n a l flow data measured in m3/s. z f c = H x t 5.1 Var (z, x t) = Var(v ) 5.2 The time subscript i s now dropped for notational convenience. It is postulated that the variance of the observation i s a function of the mean: Var(z) = a (Mean(z)) k 5.3 or log [Var(z)] = log[a] + k log [Mean(z)] 5.4 Equation 5.4 i s a linear equation with k representing the slope. It can be determined by regressing log[Var(z)] on log[Mean(z)]. Hence, the dependence of Var(z) on Mean(z) can be found. The variance and mean of the flow in eqn. 5.4 are replaced by l o c a l estimates. These estimates are moving 59 averages of order "s" centered about the current time of interest. Flow observations are made from t=1 to T. Estimates for the mean and the variance of the flow at time t are given by: s . Z z t + i Mean(z) = 1 s = z f c 5.5 2s+1 Var(z) = i=-s Z (z t + i - V* 2s 1 2s s L 3 i = -s 't + i - z 2.(2s+1) 5.6 For s=2, the averaging i s with respect to 5 data points; hence this i s termed a 5-pt. average. The smaller "s" i s , the more l o c a l the estimates are. 5.3.2 Results In a preliminary analysis, estimates are calculated with 3, 5, and 7-pt. averages. Results of the regressions show l i t t l e differences in the values obtained for k. It i s decided that a 5-pt. average w i l l be used to obtain l o c a l estimates of the mean and variance of the flow. The time period i s from A p r i l 1 to September 30. Ten years of data, 1970 to 1979 are used to determine the slope k, of eqn. 5.4. Results of the regressions are tabulated below: 60 Table 5.1 Values of the slope for years 1970-1979 year k year k 1970 2.54 1975 2.06 1971 1 .48 1976 1 .63 1972 1 .92 1977 1 .70 1 973 2.00 1 978 2.10 1974 1 .58 1979 2.64 The average of the ten k values i s 1.96. 5 .3 .3 Conclusions The averaged k value of 1.96 indicates that the variance of the measurement error i s approximately proportional to the square of the flow (eqn. 5.3). In other words, the standard deviation of the noise term i s d i r e c t l y proportional to the flow. A constant variance for the measurement error would result in k being close to zero. 5.3.4 B o x - C o x Transforms A class of transforms used for s t a b i l i z i n g the variance is known as the Box-Cox transforms. An appropriate transformation of the raw flow data can be made. The resulting series w i l l have an error term whose variance i s approximately constant with time. These tranforms are used 61 in situations where the variance can be expressed as a function of the mean M; Var(z) = 0(M) where <t> is the functional dependence. The f i r s t order Taylor Series expansion of the transformed data i s : It i s desired that Var (g(z)) be approximately constant. The functional dependence </>, determined in section 5.3.3, i s z 2 . Evaluation of eqn. 5.8 gives the appropriate transformat ion: g(z) = g(M) + g'(M)(z-M) 5.7 Var (g(z)) « g'(z) 2 0(M) 5.8 y = g(z) = ln(z) 5.9 Eqn. 5.9 indicates that the incoming flow data should be tranformed by taking the natural logarithms. The new series formed by the y's, has approximately uniform variance. 62 5.4 Maximum Likelihood Estimation of the Noise Variance 5.4.1 Theory Consider a set of T dependent observations drawn from a multivariate normal, y_ MVN(jz,a2V). The l o g - l i k e l i h o o d function, log £(y 1 r y 2 f « « . y T) = log £(y_) i s : -| log 2n - |log a2 - \ log | V | - (y_-M)'V"1 (y_-u_) 5.10 D i f f e r e n t i a t i n g the above expression with respect to a 2 and equating i t to 0 yields the maximum l i k e l i h o o d estimate of a 2. However, the evaluation of the above expression is time consuming due to | v j ; and V" 1. This is because V i s a non-diagonal matrix as the individual errors are s e r i a l l y correlated. Schweppe (1965) showed that an equivalent expression for log £(y_) can be written in terms of the innovations, v^. It i s defined as the observation forecast error, y t-y f c. These innovations are independently normally d i s t r i b u t e d with mean 0, and variance a2t^, where f f c i s calculated by the Kalman algorithm. Hence, £ i s a set of independent errors drawn from a multivariate normal, MVN(0, a 2D). D i s a diagonal matrix defined as: 63 f D = The equivalent l o g - l i k e l i h o o d expression, log £(y_) i s : Because i s a linear transformation of y^ . which i s normally d i s t r i b u t e d , equations 5.10 and 5.11 are equivalent. The forecast error , and i t s variance f^ . are quantities which can be computed by the Kalman F i l t e r at each time step, as part of the Kalman algorithm. This i s achieved by recasting the time series model into the state-space framework. Hence, the Kalman F i l t e r i s used as a tool for evaluating the l i k e l i h o o d of T dependent observations from a multivariate normal. Differen t i n g the above expression with respect to a2 y i e l d s the formula for a2. -|log 211 - | logo-2 - ^ 2 log f -t=1 T t=1 f 5.11 5.12 The quantities, and f f c are given by the Kalman F i l t e r ; where f f c = 1 + HPH'. 64 5.4.2 Results The three ARMAX models given in Chapter 4 w i l l be used for flow forecasting in Chapters 6, 7, and 8 respectively. For each model, an estimate of the noise variance based on the transformed flow data i s required. Five estimates are obtained for each model. Flow records from 1976-1980 are used, with the observations taken from A p r i l 1 to September 30. Results are given below. Table 5.2 Maximum Likelihood Estimates of the Noise Variance year Model 1 Model 2 Model AR(1)* TRANSFER ARMAX FUNCTION 1 976 .00201 .01237 .00209 1977 .00264 .01578 .00276 1 978 .00122 .01081 .00175 1979 . .00233 .01085 .00259 1980 .00287 .01888 .00303 R av .0022 .01 37 .0024 * Descriptions of these models are given in Chapter 4, section 4.2. Spe c i f i c a t i o n of the noise variance for the Kalman algorithm i s based on R_„, in subsequent chapters 6 to 8. 65 5.5 Summary The t h r e e ARMAX models d e s c r i b e d i n Chapter 4 w i l l be used f o r f l o w f o r e c a s t i n g i n C h a p t e r s 6, 7, 8 r e s p e c t i v e l y . For each model, the v a r i a n c e of the measurement e r r o r a t each time s t e p i s unknown. T h i s c h a p t e r i s concerned w i t h the e s t i m a t i o n of n o i s e v a r i a n c e . In r e a l - t i m e f o r e c a s t i n g , i t i s d e s i r e d t h a t the v a r i a n c e of the e r r o r term be t i m e - i n v a r i a n t ; a l t h o u g h i n p r a c t i c e , t h i s i s r a r e l y t r u e . E s t i m a t i o n of the n o i s e v a r i a n c e f o r each model i s d i v i d e d i n t o two p a r t s and the c o n c l u s i o n s a r e : 1. R e g r e s s i o n s of l o g [ V a r ( f l o w ) ] on l o g [ M e a n ( f l o w ) ] i n d i c a t e t h a t the v a r i a n c e i s p r o p o r t i o n a l t o the square of the f l o w . By u s i n g the Box-Cox t r a n s f o r m , i t i s det e r m i n e d t h a t the raw f l o w d a t a s h o u l d be t r a n s f o r m e d by l n ( f l o w ) . The n o i s e v a r i a n c e of the new s e r i e s i s a p p r o x i m a t e l y u n i f o r m . 2. The method of maximum l i k e l i h o o d i s used t o e s t i m a t e t h e n o i s e v a r i a n c e s of the t r a n s f o r m e d d a t a f o r each model. The l o g - l i k e l i h o o d e x p r e s s i o n i s w r i t t e n i n terms of t h e i n n o v a t i o n s and e v a l u a t e d by u s i n g the Kalman F i l t e r . Average of the maximum l i k e l i h o o d e s t i m a t e s f o r each model i s g i v e n . 66 Table 5.3 Averaged Maximum Likelihood Estimate of the Noise Variance Model 1 .0022 Model 2 .0137 Model 3 .0024 6 . A P P L I C A T I O N : A R ( 1 ) An autoregressive model of order 1, AR(1) i s considered in this chapter. It is used to predict flows at Hope one day in advance. Three ways of formulating the AR(1) model into Kalman state-space format are investigated. Hence, three schemes are used to forecast the flow at Hope. The objective of t h i s study i s to compare forecasting performances of the three state-space formulations, based on the following three c r i t e r i a : 1. RMS error of the flow forecast, 2. maximum r e l a t i v e forecast error, and 3. number of times the forecast error i s greater than 25% of the true flow. From the study, the following conclusions are drawn. The best forecasting performance i s obtained for the scheme which uses the Kalman F i l t e r to estimate the AR c o e f f i c i e n t , a. As the estimate i s given resursively, flow prediction i s made at each time step with the most recent estimate of a. On the other hand, the Kalman F i l t e r can be used to estimate the flow d i r e c t l y . The formulation which models the flow as the state, and the error term as the observation noise r e s u l t s in the worst forecasting performance. Correction to the flow estimates are not adequate as the Kalman gain matrix i s too small. Formulating the AR process as the state equation i s equivalent to not using the Kalman F i l t e r at a l l . The forecasting performance can be improved by casting the model 67 68 in state-space form and applying the f i l t e r to estimate the AR c o e f f i c i e n t recursively. The model used to describe the streamflow phenomenon at Hope i s: q t = a q t - 1 + e t 6 ' 1 where q = ln(flow) at Hope a = autoregressive c o e f f i c i e n t efc = white noise ~ N(0,a 2) 6.1 Description of the Forecasting schemes Three schemes are used to forecast flows at Hope and are described below. The streamflow phenomenon i s represented by the AR(1) model of eqn. 6.1. 6.1.1 Properties of Scheme 1 The AR(1) model i s recast into state-space framework as follows: State eqn: afc = a t _ 6.2 Obs. eqn: qfc = q 4._ 1a t + 6.3 The AR(1) process is written as the observation equation, with a i s the state variable. Hence, the Kalman F i l t e r i s used to give recursive estimates of the AR c o e f f i c i e n t . 69 Because time-invariant states are considered, there is no noise term associated with a, and thus Q = 0. The 1-step flow prediction i s given by §t+1 = ^t^t+1/t ' "here <*t+l/t ^ s based on the l a t e s t state estimate given by the f i l t e r . Under th i s scheme, the noise variance i s required as input to the Kalman algorithm. 6.1.2 Properties of Scheme 2 The state-space representation of the AR(1) model i s : State eqn. q f c = a q t - 1 + efc 6.4 Obs. eqn. q f c = q f c 6.5 In contrast to the previous scheme, the AR(1) process is written as the state equation, with qfc as the state variable. Under th i s formulation, an estimate of a is required prior to forecasting as i t is assumed known by the Kalman algorithm. An approximate least squares estimate, a L g i s used. This estimate of a i s held fixed throughout the forecasting period. Forecasts are given by the prior estimate of the state vector, Q t + 1/ t = a L S ^ t ' 6.1.3 Properties of Scheme 3 This f i n a l scheme also uses the Kalman F i l t e r to give estimates of the flow d i r e c t l y . However, the AR(1) model is " s p l i t up" in the state-space formulation: 70 State eqn: q f c = a q t - 1 6.6 Obs. eqn: qmt q f c + efc 6.7 qmt represents the measured flow. Although q f c i s the state variable, the noise term i s modelled as the measurement noise. Again, an estimate of a is required for the algorithm and a Lg i s used. Forecasts are given by the f i l t e r as Q^+1/t = aLgQ t« This d i f f e r s from Scheme 2 in the expression for q t , as w i l l be shown in section 6.5.4. An estimate for the noise variance i s also required. Table 6.1 summarizes the properties of the three schemes. 71 Table 6.1 Properties of the forecasting schemes Scheme 1 Scheme 2 Scheme 3 AR(1) = obs. eqn. AR(1) = state eqn. AR(1) " s p l i t up" estimates of a are given recursively by the Kalman F i l t e r estimate of a2 required estimate of a required prior to forecasting equivalent to no Kalman F i l t e r i n g estimate of a required prior to forecast ing estimate of o2 required 72 6.2 General Discussion For Scheme 1, the estimate of a can be updated at each time step through the use of the Kalman F i l t e r . The forecasting performance can be affected by the input s p e c i f i c a t i o n of the Kalman algorithm. For thi s formulation, these inputs are the i n i t i a l s p e c i f i c a t i o n s for the AR c o e f f i c i e n t and the noise variance; i . e . a0, P 0» R. Schemes 2 and 3 use the Kalman F i l t e r to give estimates of the flow d i r e c t l y , as q f c i s the state variable. An estimate of a i s required prior to forecasting as i t i s assumed known by the Kalman algorithm. The main disadvantage of these schemes is that a i s held fixed during the forecasting period. Hence, the performance of schemes 2 and 3 are not expected to be better than that of Scheme 1. Under Scheme 2, the 1-step flow prediction i s given by the prior estimate of the state q t + 1 y t = a L Sq f c. The lat e s t estimate of the flow i s : q t = $ t / t_. + K t[q t-q t] 6.8 The expression for the Kalman gain i s : K t = P t / t - 1 H ' t H P H ' + R l " 1 6.9 where P T / T _ 1 = a ^ P ^ + Q 6.10 Spe c i f i c a t i o n of Q aff e c t s Kfc, and thus q t . In thi s a p plication, H=1 and R=0 leads to a constant gain matrix, Kt=1. Therefore, regardless of the value for Q, Kfc i s the same even though p t / t _ i changes. Because K=1, equation 6.8 i s equivalent to: 73 6.11 Hence the 1-step forecast i s : ^t+1/t = aLS q t 6.12 This expression for the 1-step forecast could be obtained before the AR(1) model was cast into state-space form. Hence, for th i s situation (K f c=l), Scheme 2 i s equivalent to not using the Kalman F i l t e r . The s e n s i t i v i t y of t h i s scheme's performance is studied with respect to sp e c i f i c a t i o n of a. Under Scheme 3, the Kalman gain i s not equal to 1. In this case, the performance of the f i l t e r i s expected to depend on the sp e c i f i c a t i o n of q 0 , P 0, R and a. 6 . 3 Experimental Procedure The forecasting performance of each scheme i s investigated. The data used are streamflow measurements at Hope, from A p r i l 1 to September 30, for the years 1981 1982 1983. Their performance is measured by three indicators (PI): 1 . PI 2 74 This r e l a t i v e root mean square error (RRMS) expresses the average magnitude of the forecast error as a f r a c t i o n . It i s measured in terms of the actual units; i . e . y i s in m3/s. 2. PI 2 = maximum of the absolute r e l a t i v e error, 3. P I 3 = number of times the absolute flow forecast error error > 25% of the actual flow. A l l three indicators are such that the smaller their values, the better the performance. To summarize, PI, represents the average error, PI 2 denotes the maximum error, and PI 3 indicates the frequency of poor forecasts. Schemes 2 and 3 require estimates for the AR c o e f f i c i e n t prior to forecasting. The method of least squares i s used to obtain , by treating the time series model as a regression model. An approximate least squares estimate (LSE) for a is the sample c o r r e l a t i o n c o e f f i c i e n t between q t and Q t_i« Taking the dependent variable as q f c, and the independent variable as Q t _ i r five regressions are done and the sample co r r e l a t i o n c o e f f i c i e n t i s obtained in each case. Data used are streamflow records at Hope from A p r i l 1 to September 30 for the years 1976-1980. The average of the f i v e estimates i s used as a for forecasting. 75 Schemes 1 and 3 require the noise variance to be known. Si m i l a r l y , this has to be estimated prior to forecasting and the method of maximum l i k e l i h o o d i s used in thi s case. The same flow records as above are used. Assuming normality for the observations, the approximate maximum l i k e l i h o o d estimates (MLE) of a2 were obtained in Chapter 5. Again, the average of the five estimates is used as R. 6.4 Results 6.4.1 Estimate of AR c o e f f i c i e n t The LSE of a obtained for the years 1976-1980, i s presented below. Table 6.2 Least Squares Estimates of the autoreqressive c o e f f i c i e n t year 1976 1 977 1 978 1 979 1 980 average It i s noted that the estimates of a are close to 1. aLS .9963 .9928 .9959 .9969 .9945 .995 76 6.4.2 Estimate of the noise variance The MLE of a2 obtained in Chapter 5 for the years 1976-1980 are given. Table 6.3 Maximum Likelihood Estimates of the noise variance year MLE(a2) 1976 .00201 1977 .00264 1978 .00122 1979 .00233 1980 .00287 average .0022 The average value i s rounded to .002 for use in the Kalman algorithm. 6.4.3 Performance Indicators for Scheme 1 The f i l t e r ' s performance i s tested with respect to a0, P 0, and R. The range of parameter values used are given below, followed by the r e s u l t s . 77 Table 6.4 Range of values for input parameters Parameter lower l i m i t upper l i m i t a 0 .1 2 P 0 .1 10 R 2X10- 6 2 Table 6.5a Values of PI, for Scheme 1 ao=1.0 Po=3.0 R=.002 unless otherwise stated Parameter value 1981 1982 1983 a l l conditions 4% 5% 4% R=2X10"6 4% 5% 5% Table 6.5b Values of PI 2 for Scheme 1 Parameter value 1981 1982 1983 a l l conditions 17% 20% 14% R=2X10-6 18% 20% 12% 78 Table 6.5c Values of PI 3 for Scheme 1 Parameter value 1981 1982 1983 a l l conditions 0 0 0 R=2X10" 6 0 0 0 Results show that the f i l t e r i s very robust to sp e c i f i c a t i o n s of a 0 and P 0 for a l l three performance indica t o r s . It is noted that in this a pplication, the f i l t e r i s also insensitive to the s p e c i f i c a t i o n of R. 6.4.4 Performance Indicators for Scheme 2 This scheme is equivalent to forecasting without the Kalman F i l t e r . It was shown in section 6.2 that forecasts do not depend on the value of the noise variance in this a p p l i c a t i o n . The forecasting performance i s examined with respect to s p e c i f i c a t i o n of a. The range of values used are: Table 6.6 Range of values for a Parameter a lower l i m i t .95 upper l i m i t 1 .05 79 Table 6.7a Values of PI, for Scheme 2 Parameter value 1981 1982 1983 a=.95 34% 35% 34% a=.99 9% 10% 9% a=1.0 4% 5% 4% a=1.0l 10% 10% 10% a=1.05 51% 52% 51% Table 6.7b Values of PI 2 for Scheme 2 Parameter value 1981 1982 1983 a=.95 44% 48% 40% a=.99 24% 28% 20% o=1.0 18% 21% 13% O=1.01 16% 23% 18% o=1.05 66% 74% 64% 80 Table 6.7c Values of PI 3 for Scheme 2 Parameter value 1981 1982 1983 a=.95 182 179 182 a=.99 0 1 0 a=1.0 0 0 0 o=1.01 0 0 0 a=1.05 179 180 182 A l l three performance indicators show that the forecasting performance is sensitive to s p e c i f i c a t i o n of a. The best performance shows an average RRMS error (PI,) of 5%, while the worst performance (for a=1.05) results in 50%. The forecasting period spans over 183 days, from A p r i l 1 to September 30 of each year. Table 6.7c shows that for a=.95 or 1.05, the forecast error i s greater that 25% of the actual flow at almost every time step. 6 . 4 . 5 Performance Indicators for Scheme 3 The forecasting performance i s tested with respect to o-r Qo f po r and R. The range of values for the parameters are: 81 Table 6.8 Range of values for input parameters Parameter lower l i m i t upper l i m i t a .99 1.01 q 0 1 15 P 0 .3 30 R 2X10" 6 20 Results are presented below. a=1.0 qo= 7 Po = 3 R=.002 unless otherwise stated. Table 6.9a Values of PI, for Scheme 3 Parameter value 1981 1982 1 983 a=.99 91% 92% 92% a=1 .0 46% 48% 44% a=1.01 24030% 23605% 23234% R=20 46% 49% 45% R=2 46% 48% 44% R=2x10"6 46% 48% 44% a l l values of q 0 and P 0 46% 48% 44% 82 Table 6.9b Values of PI 2 for Scheme 3 Parameter value 1981 1982 1983 a=.99 1 00% 100% 1 00% a=1.0 1 30% 77% 109% a=1.01 80625% 80848% 72877% R=20 121% 78% 101% R=2 1 29% 77% 108% R=2x10"6 131% 77% 1 10% a l l values of q 0 and P 0 1 30% 77% 1 09% Table 6.9c Values of PI 3 for Scheme 3 Parameter value 1 981 1982 1983 a=.99 174 1 75 1 77 a=1 .0 1 16 1 30 1 33 a=1.01 1 78 174 1 72 R=20 1 22 1 25 1 40 R=2 1 1 5 129 1 32 R=2X10 - 6 1 16 130 1 32 a l l values of q 0 and P 0 116 130 133 A l l three Performance Indicators show that t h i s i s the worst forecasting scheme as even the smallest value for PI, is approximately 50%. Although the performance i s robust to 83 s p e c i f i c a t i o n of q 0 , P 0, and R; i t is very sensitive to the s p e c i f i c a t i o n of a. From Table 6.9a, a=1.0 results in the smallest PI,. However, using a=.99 doubles the average RMS error. Moreover, using a=1.0l leads to f i l t e r divergence. 6.5 Discussion of Results 6.5.1 Forecasting Performance of Scheme 1 For Scheme 1, the f i l t e r ' s robustness to the AR c o e f f i c i e n t , a i s not surprising. It was found in Chapter 3 that s p e c i f i c a t i o n of x 0 and P 0 has l i t t l e e f fect on the forecasting performance. It is also noted that for t h i s application, the f i l t e r i s also robust to s p e c i f i c a t i o n of the noise variance. Hence, this formulation of the Kalman F i l t e r i s most at t r a c t i v e in real-time forecasting, because a) the forecasting performance i s good; and b) the f i l t e r i s robust to i n i t i a l s p e c i f i c a t i o n of model c o e f f i c i e n t s . Indeed, the main advantage of t h i s formulation is that a can be continually adjusted by using the incoming flow data. Table 6.2 shows that the least squares estimates of a are close to 1. In fact, they are constrained to be less than 1 as they are the sample co r r e l a t i o n c o e f f i c i e n t s . However, the Kalman F i l t e r does not have th i s l i m i t a t i o n . Consequently, the Kalman F i l t e r can be used to estimate the AR c o e f f i c i e n t for non-stationary systems where a i s greater than 1. For t h i s application, the estimates of a given by 84 the Kalman F i l t e r are also close to 1. Hence, they support the regression estimates. 6.5.2 Forecasting Performance of Scheme 2 The AR(1) model i s written as the state equation of the Kalman model. For K=1, t h i s i s equivalent to forecasting without the Kalman F i l t e r . An estimate of a i s required prior to forecasting. Results show the f i l t e r ' s s e n s i t i v i t y to s p e c i f i c a t i o n of the AR c o e f f i c i e n t . Values of a outside a 5% range with respect to a L g results in larger PI values. In practice, a cannot be estimated with an accuracy of a few percent. Even s l i g h t seasonal v a r i a t i o n can cause the true value of a to wander by t h i s amount. 6.5.3 Forecasting Performance of Scheme 3 Forecasting under Scheme 3 r e s u l t s in the worst performance. Large increases in a l l Performance Indicators are noted when s p e c i f i c a t i o n of a i s changed by as l i t t l e as 1% from a L g . Therefore, t h i s i s not a p r a c t i c a l forecasting procedure, as a cannot be estimated with t h i s degree of accuracy. Poor estimates of a once chosen cannot be corrected as i t i s held fixed during the forecasting period. Flow predictions are given by: $t+l/t = aLS $t 6 * 1 3 q. i s the Kalman estimate of the state vector for time t: $t = V t - 1 + K t [ ( 3 m t " q i \ ] 6.14 85 or $t = q t / t - i ^ - V + K t q m t 6 - 1 5 There are two l i m i t i n g cases of in t e r e s t : 1. K t = 1 results in flow forecasts given by q t + 1 y t = a L Sqm t. This i s equivalent to not using the Kalman F i l t e r (Scheme 2). 2. Kfc=0 results in flow forecasts given by Q t + i / t = a L s q t / t - 1 * T n i s corresponds to not using the measurement information at time t. This scheme leads to propagation of errors as t increases. In fact, t h i s i s the reason why t h i s forecasting scheme yield s such poor re s u l t s . Successive flow predictions are based on the previous prediction without consideration of the measurements. This then explains why a=1.0l results in much larger Performance Indicators than a=.99. Flow forecasts in this situation is approximated by q 0 raised to an exponent a f c, where q 0 can be taken to be the i n i t i a l flow. For large t, 1.01*" r e s u l t s in forecasts, which deviate from the actual flow by a greater amount than .99*". It i s shown below that for t h i s a p p l i c a t i o n , Kfc i s approximately 0. For H=1, the expression for the Kalman gain reduces to: K t " P t / t - 1 / [ P t / t - 1 + R ] 6 ' 1 6 K t i s small i f p t / t _ i i s much less than R. Starting with P 0, Pyo = a 2 L s P ° + °- ~ p ° • 86 Hence, K, ~ p 0/[ p 0+R], which gives P, ^ ( I - K ^ P Q . 1. If a large P 0 i s spe c i f i e d , then K,~1 which results in P^O. Subsequent P F C w i l l a l l be close to 0. 2. If P 0 i s small compared to R, then K^O. However, P,~P 0 which i s s t i l l small. Hence after a couple of time steps, p t / / t _ i approaches 0 regardless of P 0 . This situation should be distinguished from Scheme 1 where the AR process is written as the observation equation. In that case, a2 = R also. The state variable however, i s the AR c o e f f i c i e n t . Hence the Kalman gain i s applied to a. Since small adjustments to a are adequate for improving the observation forecast, only a small Kalman gain i s necessary. If the state variable i s the flow (as in Scheme 3) then small Kfc values are not s u f f i c i e n t to correct the flow estimates. Therefore, this forecasting scheme i s the worst one considered in t h i s chapter because: 1. No adjustment can be made to a. 2. Flow prediction at each time step i s not corrected adequately as K i s too small. 87 6 . 6 Summary The AR(1) model is used to describe streamflow phenomenon at Hope. Results for the performance of the forecasting schemes are given below. Scheme 1 i s the best forecasting procedure. This corresponds to formulating the AR(1) model in the observation equation with a as the state v a r i a b l e . Hence, thi s allows continual updating of the AR c o e f f i c i e n t . The scheme i s best both in terms of flow forecasts obtained, and robustness to i n i t i a l s p e c i f i c a t i o n of a. Scheme 2 formulates the AR(1) process as the state equation in the Kalman model. This is equivalent to not using the f i l t e r at a l l . Although the performance can be as good as Scheme 2, forecasts are sensitive to the s p e c i f i c a t i o n of the AR c o e f f i c i e n t . Scheme 3 i s the worst forecasting procedure. Not only are the flow predictions sensitive to a, the performance i s always poor. Under th i s scheme, a small Kalman gain leads to i n s u f f i c i e n t correction to qJ. + 1/4.. 7. APPLICATION: TRANSFER FUNCTION MODEL The AR(1) model in chapter 6 gives the 1-step ahead forecast for the flow at Hope. Application of the Kalman F i l t e r does not necessarily y i e l d improved flow predictions. The best forecasting scheme i s to use the Kalman F i l t e r to estimate the AR c o e f f i c i e n t recursively. This allows flow predictions to be made with the la t e s t estimate, a f c. The formulation corresponds to writing the AR(1) model in the observation equation. Often in flood management however, forecasts are required for several days in advance. There i s a problem with using the AR(1) model in a Kalman F i l t e r framework to predict observations more than 1 step ahead. In the Kalman model, Hfc which relates the observation to the state, i s considered known at time t. For the 1-step forecast, Ht + 1 = Q t« But for the 2-step forecast, H f c + 2 = q t + 1 « The problem is that tomorrow's flow is unknown. In practice, the 1-step forecast i s used as an estimate. However, the variance of this forecast error (<3t +2 ~ Q t + 2^' should not be calculated using the Kalman algorithm. One alternative i s to use a di f f e r e n t s t a t i s t i c a l model, which avoids t h i s problem, to predict the flows 2 days in advance. A transfer function model i s used to relate flows from upstream stations to the flow at Hope. This contrasts with the time series model of Chapter 6. In this case, the independent variables do not involve past values of the flow at Hope, the dependent vari a b l e . This input-output model 88 gives the 2-day ahead forecast at d a i l y i n t e r v a l s , as measurements are received everyday. The model used to describe the flow at Hope i s a transfer function model with constant c o e f f i c i e n t s : q H t = [ q T t _ 2 q N t _ 2 ] I + e t 7. where q = ln(flow) §_ = vector of c o e f f i c i e n t s e = white noise~N(0,a 2=R) station superscripts: H = Hope T = Texas Creek N = Near Spences Bridge Stations T and N are located upstream of Hope. Their geographic locations are shown in Chapter 1 , F i g . 1 . 1 . The objective i s to determine whether flow forecasts given by th i s model can be improved by using the Kalman F i l t e r to estimate §_ in real-time. Performance of the forecasting schemes are based on the same indicators used for the A R ( 1 ) study. The following conclusions are drawn: 1. Flow forecasts are not improved s i g n i f i c a n t l y by the Kalman F i l t e r to be of p r a c t i c a l i n t e r e s t . 90 2. However, under the Kalman F i l t e r scheme, forecasting performance i s robust to s p e c i f i c a t i o n of £. 7.1 Description of Forecasting Schemes Only two forecasting schemes are considered in this study. In Chapter 6, three state-space formulations of the AR model were considered. It was found that the formulation which models the AR c o e f f i c i e n t as the state y i e l d s improved forecasts. Hence, this i s the only state-space formulation considered in th i s study. 7.1.1 Properties of Scheme 1 No Kalman f i l t e r i n g i s employed in the f i r s t scheme. The transfer function model assumes that the c o e f f i c i e n t s are time-invariant. Flow 2 days in advance i s given by the 1-step forecast: q H t + 2 - [ Q T t q N t l I L S 7.2 where j 3 L S i s the least squares estimate obtained prior to forecasting. 7.1.2 Properties of Scheme 2 Kalman F i l t e r i n g i s used and the transfer function model i s cast into state-space format as follows: 91 = "»'t-r = _»'t-i. 7.3 q H t = [ Q T t _ 2 q N t _ 2 ] lt + e f c 7.4 The c o e f f i c i e n t s j31, |32, are modelled as time-invariant states. Therefore, the Kalman F i l t e r i s used to give estimates of £ f c recursively at each time step. The flow 2 days in advance i s : q " t + 2 - [ q T t q N t ] lt 7.5 where |?t i s based on the l a t e s t estimate given by the Kalman F i l t e r . In order to apply the Kalman algorithm, an estimate of the noise variance is required. 7.2 General Discussion Results of the previous study imply that the scheme which does use the Kalman F i l t e r would y i e l d better flow forecasts. This i s because £ i s estimated recursively in real-time. Consequently, flow predictions are made at each time step with the latest estimate of the c o e f f i c i e n t s . In addition, i t i s expected that under Scheme 2, the forecasting performance would be robust to the i n i t i a l s p e c i f i c a t i o n , j 3 0 . This i s not l i k e l y to be true for Scheme 1. Thus, performance of Scheme 1 i s examined with respect to d i f f e r e n t values of j 3 r c : « Scheme 2 i s examined with respect 92 to input s p e c i f i c a t i o n of J30 and R. It i s not known how the forecasting performance of t h i s proposed model compares with that of the autoregressive model. This i s a question of model i d e n t i f i c a t i o n and i s not the main objective of t h i s study or thesis. However, the following reasoning can be given. Stations with large drainage areas such as Hope, have r e l a t i v e l y slow system response. This is due to the storage c h a r a c t e r i s t i c s of large basins. Quimpo (1973) showed that t h i s phenomenon can be represented by an autoregressive process. With t h i s consideration, i t would not be surprising to find that the forecasting performance of the AR(1) model i s better than that of the transfer function model. 7.3 Experimental Procedure The forecasting performance of each scheme is investigated using streamflow measurements from three stations: Hope, Texas Creek, and Near Spences Bridge. The forecasting period i s from A p r i l 1 to September 30 for 1981 1982 and 1983. The performance of each scheme i s measured by the three indicators as defined in Chapter 6. These are: 1. PI, which represents the average error 2. PI 2 which denotes the maximum error 3. PI 3 which indicates the frequency of poor forecasts Small PI values are i n d i c a t i v e of good forecasting performance. 93 Flow p r e d i c t i o n under Scheme 1 r e q u i r e s an e s t i m a t e f o r J3 p r i o r t o f o r e c a s t i n g . The method of l e a s t squares i s used t o o b t a i n £ Lg. Flow r e c o r d s used f o r t h i s e s t i m a t i o n a r e based on f i v e y e a r s of d a t a , 1976-1980. The same time p e r i o d of the year i s used; A p r i l 1 t o September 30. An average of the f i v e l e a s t squares e s t i m a t e s (LSE) i s used as J3 f o r f o r e c a s t i n g . R e a l - t i m e f o r e c a s t i n g w i t h the Kalman F i l t e r r e q u i r e s an e s t i m a t e f o r the n o i s e v a r i a n c e , R. Approximate maximum l i k e l i h o o d e s t i m a t e s (MLE) were o b t a i n e d i n Chapter 5 f o r the y e a r s 1976-1980. A g a i n , the average of the s e f i v e e s t i m a t e s i s used as s p e c i f i c a t i o n f o r R. 7.4 R e s u l t s 7.4.1 E s t i m a t e of the model c o e f f i c i e n t s R e g r e s s i o n a n a l y s i s g i v e s v a l u e s of p 2 c l o s e t o 0 f o r a l l y e a r s , 1976-1980. Hence o n l y the l e a s t s q uares e s t i m a t e s f o r |3, a r e g i v e n below. 94 Table 7.1 Least Squares Estimates of the model c o e f f i c i e n t year ^LS 1976 (1 .0572) 1977 (1.0514) 1978 (1.0745) 1 979 (1.0556) 1980 (1.0675) average (1.0613) The regression results indicate that flows from station N do not have a s i g n i f i c a n t contribution in predicting the 2-day advance flow at Hope. Because j32 = 0, only the station at Texas Creek is used as the upstream input in thi s input-output model. Thus, p ° A V = 1.0613. 7 . 4.2 Estimate o f the noise variance The f i v e MLE of a2 from Chapter 5 are given below. 95 Table 7.2 Maximum Likelihood Estimates of the noise variance year R 1976 .01237 1977 .01578 1978 .01081 1979 .01085 1980 .01888 average .0137 R i s rounded off to .015 for the Kalman algorithm input 7.4.3 Performance Indicators for Scheme 1 Results of the performance measures using j3 A V are given below. Table 7.3 Performance Indicators for Scheme 1 PI i PI 2 PI 3 1981 10% 31% 3 1982 1 7% 48% 24 1 983 18% 41% 35 The analagous scheme of the AR(1) model i s compared. 96 This corresponds to flow prediction using a L g without the Kalman F i l t e r . Table 7.4 Performance Indicators for the AR(1) model 1981 1982 1983 PI , 5% 7% 6% PI 2 21% 24% 16% PI 3 0 0 0 Values of the PI for the transfer function model (Table 7.3) are higher than those of the AR(1) model (Table 7.4). This indicates a l e v e l a performance which is not as good as that obtained with the autoregressive model. The s e n s i t i v i t y of Scheme 1 with respect to input s p e c i f i c a t i o n i s examined. Streamflow data from 1982 are used as t h i s resulted in the worst performance. Table 7.5 S e n s i t i v i t y of the Performance Indicators to the model c o e f f i c i e n t 0=1.0613 0=1.04 PI, 17% 20% PI 2 48% 46% PI 3 24 37 97 The f i r s t value of 0 i s 0 Ay. The second value of 0 i s an ar b r i t r a r y choice which represents a 2% change from the averaged f3A V» There i s an increase in PI, (average error) and PI 3 (frequency of poor forecasts) for 0=1.04. 7 . 4 . 4 Performance Indicators for Scheme 2 » Values of the performance indicators for Scheme 2 are presented: Table 7.6 Performance Indicators for Scheme 2 P I 1 P I 2 P I , 1981 9% 40% 3 1 982 1 9% 54% 32 1983 1 3% 33% 12 Tables 7.3 and 7.6 are compared. For the 1981 data, comparable forecasting performance i s obtained for the two schemes. In both cases, the number of times the forecast error i s in excess of 25% of the actual flow, (PI 3) i s minimal. The 1982 data reveal a s l i g h t l y larger difference in the performance indicators between the two schemes. A l l three performance indicators are larger for Scheme 2 which uses the Kalman F i l t e r . This implies that better flow forecasts can be obtained without continual updating of £. 98 However, performance indicators for the 1983 data indicate that better forecasts are obtained i f the Kalman F i l t e r i s used. For t h i s scheme, the s e n s i t i v i t y of i t s performance with respect to the i n i t i a l s p e c i f i c a t i o n of j3, and R i s examined using the 1982 data. The range of values used for these parameters are: Table 7.7 Range of Values for input parameters Parameter value lower l i m i t upper l i m i t i n i t i a l £ (-.1.9) (1.1) R .0015 15 Under a l l input s p e c i f i c a t i o n s , the PI values are the same. Table 7.8 Performance Indicators for Scheme 2 PI 1 PI 2 PI 3 1 9% 54% 32 This contrasts from Table 7.5 where the PI values under Scheme 1 are sensitive to the s p e c i f i c a t i o n of j3. 99 7.5 Discussion of Results Performance Indicators for the two schemes do not c l e a r l y indicate which scheme gives better forecasts. However, results of the s e n s i t i v i t y analyses favor the scheme which uses the Kalman F i l t e r . This i s because the forecasting performance under Scheme 2 i s robust to i n i t i a l s p e c i f i c a t i o n of the model c o e f f i c e n t s . It i s found that without the Kalman F i l t e r , a marked decrease in performance can occur. This was i l l u s t r a t e d with a 2% change in the s p e c i f i c a t i o n of 0 for the 1982 dataset. A regression performed on the 1982 data found that | 3 L S = 1 .062. This i s very close to 0 A V of 1.0623. Hence, results in Table 7.3 represent the best performance possible for the 1982 data under Scheme 1. 100 7 . 6 Summary A time-invariant transfer function model i s used to describe the streamflow phenomenon at Hope. The independent variables are flows from two upstream stations lagged 2 days behind that of Hope. It was found that: 1. The forecasting performance of t h i s model i s not s i g n i f i c a n t l y improved or degraded by applying the Kalman F i l t e r to recursively estimate J3. 2. However, i t i s recommended that the Kalman F i l t e r be used because the forecasting performance i s robust to sp e c i f i c a t i o n of This i s important in situations where the model c o e f f i c i e n t s cannot be estimated accurately or i f they change slowly through time. 8. APPLICATION: ARMAX MODEL The Kalman F i l t e r has been applied to autoregressive and transfer function models to y i e l d estimates of the model c o e f f i c i e n t s recursively. These models are used for predicting flows at Hope 1 and 2 days in advance respectively. The advantage of using the Kalman F i l t e r i s that the forecasting performance i s robust to i n i t i a l s p e c i f i c a t i o n of the model c o e f f i c i e n t s . As these are usually unknown in hydrologic applications, t h i s feature i s appreciated by practicing engineers. In addition, use of the Kalman F i l t e r gives improved flow forecasts for the AR ( 1 ) model. An ARMAX model for describing the flow at Hope is considered in this chapter: q H t - c i q H t _ l + C z q T t - 2 + c 3 S N T - 2 + e t 8 ' 1 This i s an AR ( 1 ) with upstream inputs. The autoregressive part of t h i s model represents storage c h a r a c t e r i s t i c s of the river basin, while the exogeneous inputs represent contributions from upstream stations. In this study the Kalman F i l t e r i s applied to the ARMAX model to give estimates of the ARMAX c o e f f i c i e n t s resursively. Flow predictions 1 and 2 days in advance are made. The performance i s therefore, measured for both the 1-step and 2-step flow forecasts. The objective i s to compare the 1-step forecasting performance of the ARMAX model to that of 1 0 1 1 0 2 the AR ( 1 ) model and the 2-step performance to that of the regression model. Results of the study indicate that the 1-step forecasting performance of the ARMAX model i s comparable to that of the AR ( 1 ) model. In addition, the 2-step performance is better than that of the transfer-function model. 8.1 Properties of the Forecasting Scheme The state-space formulation of the ARMAX model i s : C i = "c r c 2 = c 2 c 3 = c 3 H r H T N -. , q t - [q t_! q t _ 2 q t - 2 ] £ t + e t The Kalman F i l t e r i s used to give estimates of the model c o e f f i c i e n t s c,, c 2 , and c 3 at each time step. Performance is measured for both the 1-step and 2-step flow forecasts. The 1-step forecast is given by: q H t + 1 = [ q H t q T t _ ! q N t _ i ] £ t+l/t 8 The 2-step forecast is given by: q t + 2 = fq t + 1 q t q t J S.t+2/t 8'5 where £ t + 1 y t and £ t + 2 / t a r e ^ a s e < ^ o n t n e l a t e s t state 103 estimate given by the f i l t e r . As q i s unknown, i t i s replaced by the 1-step forecast, q t + 1 • The autoregressive nature of this model can result in the problem of unknown system matrix, H. This occurs when predicting flows more than one time step ahead. 8.2 General Discussion It i s desired in practice that the 1 and 2-step performance of the ARMAX model be better than or equal to that of the AR(1) and transfer function models respectively. F i r s t of a l l , only one model needs to be used. Secondly, there i s only one noise variance to estimate. The station at Hope i s characterized by a large drainage area. Thus, i t i s expected that the autoregressive part of the combined model w i l l dominate. This means that estimates of c, w i l l be r e l a t i v e l y larger than those of c 2 and c 3 . With this reasoning, the 1-step forecasting performance should be similar to that of the AR(1) model. The Kalman algorithm requires s p e c i f i c a t i o n of c 0 , P 0, and R. From the studies of Chapters 6 and 7, i t has been concluded that forecasting performance i s robust to i n i t i a l s p e c i f i c a t i o n of the model c o e f f i c i e n t s . Therefore, the s e n s i t i v i t y of the f i l t e r ' s performance i s examined with respect to the s p e c i f i c a t i o n of the noise variance only. One issue of concern here i s that the 2-step prediction involves unknown elements in the t r a n s i t i o n matrix, Hfc+2 of the Kalman model. For the ARMAX model, Hfc+2 = 104 H T N fQ t + 1 Q t Q tl« T n e Kalman algorithm assumes that t h i s i s known. As tomorrow's flow at Hope is unknown at time t, i t s 1-step forecast i s used as an estimate. V i o l a t i n g t h i s assumption results in an increase in the variance of the forecast error. Therefore, the 2-step forecasting performance of the f i l t e r gives an indication of the robustness of the f i l t e r to the problem of unknown H in a pr a c t i c a l s i t u a t i o n . 8 . 3 Experimental Procedure The forecasting performance i s measured by the three indicators described in Chapter 6. For thi s ARMAX model, the performance i s measured with respect to both the 1-step and 2-step forecast errors: 1. 1-step forecast error: Q t + 1/ t _ P^+i 2. 2-step forecast error: ^t+2/t ~ qt+2 The indicators are based on these forecast errors measured in actual units of m3/s. These are: 1. PI, which represents the average error 2. PI 2 which denotes the maximum error 3. PI 3 which indicates the frequency of poor forecast. The nature of these indicators i s such that the smaller the values, the better the performance. The forecasting period is from A p r i l 1 to September 30 for 1981 1982 1983. 105 For the Kalman algorithm input, an estimate of the noise variance i s required, and the approximate maximum li k e l i h o o d estimates obtained in Chapter 5 are used. 8.4 Results 8.4.1 Estimate of the noise variance Approximate maximum li k e l i h o o d estimates (MLE) of a2 were obtained in Chapter 5 for the years 1976-1980. The average of these i s used as the s p e c i f i c a t i o n for R in the Kalman algorithm. Table 8.1 Maximum Likelihood Estimates of the noise variance year R 1976 .00209 1977 .00276 1978 .00175 1979 .00259 1980 .00303 average .0024 R i s rounded to .0025 for the Kalman s p e c i f i c a t i o n . 106 8.4.2 Performance indicators for the ARMAX model Two sets of performance indicators (PI) are given. The f i r s t set refers to the 1-step performance of the Kalman F i l t e r . The second set refers to the 2-step performance. 1-step Performance Indicators Table 8.2a Variation of PI, to noise s p e c i f i c a t i o n Value of R 1981 1982 1983 .025 4% 4.5% 3.5% .0025 4% 4.5% 3.5% .00025 4% 4.5% 4.0% Table 8.2b Vari a t i o n of PI 2 to noise s p e c i f i c a t i o n Value of R 1981 1982 1983 .025 26% 17% 21% .0025 34% 17% 24% .00025 35% 17% 24% 107 Table 8.2c Variation of PI 3 to noise s p e c i f i c a t i o n Value of R 1981 1982 1983 .025 1 0 0 .0025 1 0 0 .00025 1 0 0 2-step Performance Indicators Table 8.3a Variation of PI, to noise s p e c i f i c a t i o n Value of R 1981 .025 9% .0025 11% .00025 11% 1 982 8.5% 8.5% 8.5% 1 983 7% 7% 7.5% Table 8.3b Variation of PI 2 to noise s p e c i f i c a t i o n Value of R 1981 .025 62% .0025 104% .00025 108% 1 982 32% 32% 32% 1 983 33% 34% 44% 108 Table 8.3c Variation of P I 3 to noise s p e c i f i c a t i o n Value of R 1981 1982 1983 .025 2 2 .0025 2 2 .00025 3 2 Tables 8.2a and 8.3a show that values of PI, are less than 12% for both the 1 and 2-step performance. They also indicate that the r e l a t i v e RMS error for the 2-step forecasts i s twice that of the 1-step forecasts. Tables 8.2(a,b,c) show the r e l a t i v e i n s e n s i t i v i t y of the 1-step PI to s p e c i f i c a t i o n of R. Table 8.3b indicates a larger decrease in P I 2 for the 2-step forecasts when R i s specified larger. 8.4.3 Comparison of 1-step Performance Indicators The 1-step forecasting performance of the ARMAX model is compared to that of the AR(1) model. In both cases, the Kalman F i l t e r i s used for updating the model c o e f f i c i e n t s . 1 09 Table 8.4 Values of PI for ARMAX and AR(1) models PI:ARMAX 1981 1982 1983 :AR(1) PI , 4% 5% 4% 4% 5% 4% PI 2 34% 17% 24% 17% 20% 14% PI 3 1 0 0 0 0 0 Values of PI, and PI 3 are approximately the same under both models. However, the ARMAX model yiel d s larger values for PI2# i . e . larger maximum errors are obtained. 8.4.4 Comparison of 2-step Performance Indicators The 2-step forecasting performance i s compared to that of the transfer function model. The formulation used corresponds to l e t t i n g the model c o e f f i c i e n t s be the state vector in the Kalman model. 110 Table 8.5 Values of PI for ARMAX and Transfer Function models PI:ARMAX 1981 1982 1983 :TRANSFER FUNCTION PI, 11% 8.5% 7% 9% 19% 13% PI 2 104% 32% 34% 40% 54% 33% PI 3 2 1 2 3 32 12 Values of PI, for the ARMAX model are comparable or less than those for the input-output model. Results for PI 2 are more varied. Forecasting under the ARMAX model yiel d s a higher P I 2 value for the 1981 data; while a lower PI 2 value i s obtained for 1982 data. The 1983 data results in almost the same P I 2 values for both models. For P I 3 , the ARMAX model y i e l d s lower values in a l l three cases. In fact, the number of times the 2-step forecast error i s greater than 25% of the actual flow i s at most 2. This compares with 32 for the transfer function model. Hence, the ARMAX model seldom gives poor forecasts. 111 8.5 D i s c u s s i o n of R e s u l t s 8.5.1 Performance of the ARMAX model Tables 8.2c and 8.3c indicate that t h i s model gives r e l a t i v e l y few forecasts which result in forecast error greater than 25% of the actual flow. This i s true for both the 1-step and 2-step flow predictions. The 1-step forecasting performance i s r e l a t i v e l y i n s e n s i t i v e to sp e c i f i c a t i o n of the noise variance. However, s p e c i f i f y i n g R larger results in better 2-step performance indicators. This i s not surprising as the true noise variance for the 2-step forecast i s expected to be larger. 8.5.2 Comparison of 1-step performance The 1-step forecasts of the ARMAX model are compared with those of the AR(1) model. Results do not give a strong indication as to which model is better, as comparable values for a l l three PI are obtained. If one i s only interested in predicting the flow one day in advance, then the AR(1) model i s adequate. In fact, i t i s preferred over the ARMAX model because i t gives a parsimonious representation. 8.5.3 Comparison of 2-step performance As noted in section 8.4.4, the frequency of poor forecasts (PI 3) i s less for the ARMAX model than the transfer function model. PI 3 i s an important indicator i f ove r a l l consistency and r e l i a b i l i t y of the forecasts are 112 c r i t e r i a in practice. In addition, the average forecast error, represented by PI,, for the ARMAX model i s less than or equal to that of the input-output model. Hence, in terms of model i d e n t i f i c a t i o n , there i s a sl i g h t preference for the ARMAX model for predicting streamflow 2 days in advance. The fact that reasonable performance i s obtained for the 2-step forecasts with the ARMAX model leads to the following conclusion. For the Kalman model, ca l c u l a t i o n of the 2-step forecast requires that q i s known. As th i s i s unknown in practice, i t s estimate given by the Kalman F i l t e r as the 1-step forecast i s used. Although t h i s v i o l a t e s the assumption of the Kalman model, i t i s found in t h i s application that the 2-step forecasting performance of the Kalman F i l t e r is acceptable for engineering purposes. In fact, the performance i s better than that of the transfer function model. For the hydrologist, i t i s important that the flow predictions are given with their associated standard error. In the case of the 1-step forecast, t h i s variance i s given by the Kalman F i l t e r as HPH' + R. This assumes that H i s a deterministic known quantity. When elements of the tr a n s i t i o n matrix H are unknown, an estimate i s used; as i l l u s t r a t e d in the 2-step forecast above. Results of the 2-step PI in this study show that reasonable forecasting performance can s t i l l be obtained. However, there s t i l l remains the problem of what expression i s to be used for the variance of the 2-step forecast error. This i s addressed in 1 13 the next c h a p t e r . 8.6 Summary I f f low p r e d i c t i o n s more than 1-day i n advance are r e q u i r e d , the ARMAX can be used . T h i s i s because i t has a comparable performance as the AR(1) model . In a d d i t i o n , the 2 - s t ep f o r e c a s t i n g performance i s b e t t e r than t h a t of the i n p u t - o u t p u t model . In h y d r o l o g i c a p p l i c a t i o n s , the assumpt ion of known system m a t r i x , H i s o f t e n not s a t i s f i e d when p r e d i c t i n g f u t u r e o b s e r v a t i o n s more than one time s t e p ahead . T h i s s tudy i n d i c a t e s that u s i n g an e s t i m a t e , H s t i l l r e s u l t s in a reasonab le f o r e c a s t i n g per formance . M o r e o v e r , i t i s found tha t in t h i s p a r t i c u l a r example, the 1 and 2 - s t e p f o r e c a s t s are q u i t e i n s e n s i t i v e to the s p e c i f i c a t i o n of R. 9. VARIANCE OF THE FORECAST ERROR In flood management, the engineer i s often required to predict streamflows several days ahead. The r e l i a b i l i t y of these predictions is r e f l e c t e d by t h e i r associated standard errors. For the Kalman F i l t e r state-space model, * t = * * t - l + * t 9 * 1 Y_t = H tx t + v t 9.2 the 1-step forecast error for y_ , known as the innovation, is given. Its variance, calculated by the Kalman algorithm assumes that the system matrix H, i s known. In t h i s chapter, a general expression for the variance of the observation forecast error when both H and x are unknown, i s developed. ARMAX models have been used for streamflow modelling in th i s thesis. The state-space formulation considered here i s where the ARMAX c o e f f i c i e n t s are the state variables. The system matrix, H contains past streamflows. The k-step flow forecast, £ t + k i s given by Ht+k-t+k" F o r t h e A R ( 1 ) model, flow predictions more than 1 time step ahead requires knowledge of future flows, hence H t +^ i s unknown. For instance, ^ t + 2 = fit + 2x-t + 2 = ^-t+]-t + 2 9 , 3 In practice, the 1-step prediction i s used as an estimate for Y_t+1 . Because £ t + 1 and x.t +2 are both based on past values of they are not necessarily independent. The objective of t h i s study i s to develop a general expression for the variance of the forecast error when both H and x are unknown, and their estimates are correlated with each other. 1 14 1 1 5 From here on, the time s u b s c r i p t s a r e dropped f o r n o t a t i o n a l c o n v e n i e n c e . 9.1 Background H a r r i s o n & Stevens (1976) s u g g e s t e d t h a t the mean square e r r o r (MSE) of the k-s t e p f o r e c a s t can be o b t a i n e d from the Kalman F i l t e r . For k=1, the f i l t e r g i v e s Cov (y_-£) = HPH' + R 9.4 where R and P a r e c o v a r i a n c e m a t r i c e s f o r the o b s e r v a t i o n n o i s e , and s t a t e e s t i m a t e s r e s p e c t i v e l y . As p o i n t e d out by P r i e s t l e y i n the d i s c u s s i o n of the above paper (1976), t h i s i s not a p p r o p r i a t e i f H i s unknown. F e l d s t e i n (1971) d e v e l o p e d a fo r m u l a f o r Cov (y-y) when bot h H and x a r e unknown. R e g r e s s i o n models a re c o n s i d e r e d i n t h a t c o n t e x t . H c o r r e s p o n d s t o a d e s i g n m a t r i x which c o n t a i n s independent v a r i a b l e s , and x i s a v e c t o r of r e g r e s s i o n c o e f f i c i e n t s . The fo r m u l a g i v e n by F e l d s t e i n assumes t h a t H and x a r e independent. T h i s assumption p r e s e n t s some l i m i t a t i o n i n h y d r o l o g i c a p p l i c a t i o n s . As not e d f o r the AR(1) example, H and x a r e l i k e l y t o be c o r r e l a t e d as they a r e both based on pas t v a l u e s of y. 9.2 Problem d e f i n i t i o n The s t a t i s t i c a l model c o n s i d e r e d i s : y_ = H x + £ 9.5 1 16 where 1 vector of observations x vector of unknown parameters H system matrix e noise term The objective is to fi n d Cov (y_-£) when H and x are both unknown. Assumptions are: 1. H, x are unbiased estimators for H and x respectively. 2. The following covariance submatrices are assumed known. Z,!= E (fi-h)(fi-h)' E,2= E (fi-h)(x-x)' L22= E (x-x)(x-x)' where fi = vec H, h = vec H. The operator "vec" on matrix H stacks i t s columns one under the other, resulting in a column vector, h. The forecast £ is given by Hx. No assumption i s made on the re l a t i o n between H and x. The covariance matrix for the forecast error i s : Evaluation of Cov (Hx) involves c a l c u l a t i n g the expected value of squared quantities of H and x. This requires the fourth moment of their joint d i s t r i b u t i o n . Expression for the fourth moment can be obtained in terms of lower order moments, i f H and x are assumed to be j o i n t l y d i s t r i b u t e d as Cov (y_-£) = Cov (Hx) + Cov U) 9.6 117 a multivariate normal. A convenient way of writing products of the elements of H and x i s through the use of the Kronecker product. A preliminary result i s given in the next section before the derivation of the general formula. The resu l t contains the expressions for the covariances of products of normal random variables. It i s written in terms of the Kronecker product. 9.3 Covariances of Products of normal random variables The Kronecker product i s defined as follows: Given 2 matrices A , B ., then A ® B i s a (ms,nt) matrix mxn sxt with submatrices consisting of a^jB. Therefore, the Kronecker product of a n-dimensional vector, i s a super long vector of length n 2. From Magnus and Neudecker (1979), the following result is used: If a i s distributed as N n(y, V) where "n" i s the length of a, then Cov (a ® a) i s given by (I + K n) (V ® V + V ® MM' + MM' ® V) 9.7 K i s a n 2 by n 2 matrix defined such that K vec(A) = n J vec(A'). For the purpose of the thesis, the following derivation i s r e s t r i c t e d to a scalar model for y f c. The AR(1) model i s considered where y f c = Y t_ 1 + e^ .. Y^-i a n o - a t c o r r e s P o n c l to Hfc and xfc respectively. Hence, both H and x are scalars. Flow forecasts greater than 2 steps ahead, require estimates 1 18 for H. The expression for the variance of the forecast error i s developed in the next section. 9.4 Variance of the forecast error for the scalar case Identifying the vector a from section 9.3 as the Kronecker product, a ® a = fl ft H x x fi X X 9.9 The scalar quantity Var (fix) i s given by the 2nd (or 3rd) element of the above vector product. It i s assumed that fi N , Z , , L , 2 Z 1 2 2-2 2 where Z^j have been defined in section 9.2. In this development, they are scalar quantities. Result 9.7 from section 9.3 i s evaluated below, K , = H "2, 1 Z ! 2 i d e n t i f i e d as X 2 Z 2 2 1 0 0 0 2 0 0 0 0 0 1 0 (I + K 2) = 0 1 1 0 0 1 0 0 0 1 1 0 0 0 0 1 0 0 0 2 119 The following three matrices are symmetric r 2 ^ 1 i Also V ® V = 1^12 1^12 Z 1 2 1^22 z 2 1 2 Z ! 2 Z 1 Z 2 2 Z ! 2 Z Z 2 2 2 2 2 2 2 V ® M M ' = H ^ 1 1 HxZ, 1 H 2Z 1 2 HxZ j 2 x 2Z, 1 HxZ 1 2 X 2Z ! 2 H 2Z 2 2 HxZ 2 2 X 2 Z 2 2 H L 1 1 H2Z, 2 HxZ 1 1 HxZ ! 2 H 2Z 2 2 HxZ 1 2 HxZ 2 2 M M ' ft V = x 2 Z 1 i x 2 Z 1 2 X 2 Z 2 2 Addition of these 3 terms and premultiplying by (I + K 2) gives the expression for Var(Hx): H 2Z 2 2 + Z „ (x 2 + Z 2 2) + £ 1 2 < 2 H x + ^ 2 ) 9.9 Since, Var (y-y) = Var e + Var (Hx) the variance of the forecast error i s : (R + H 2Z 2 2) + Z V 1 (x 2 + Z 2 2 ) + Z 1 2 (2Hx + Z 1 2 ) 9.10 1 20 Eqn. 9.10 gives the general formula for the variance of the forecast error when both H and x are unknown scalar quantities. Special cases can be obtained and are discussed below. 1. If H i s unknown, but i t s estimate H is uncorrelated with x, then Z 1 2 = 0. This assumption is v a l i d for regression models where the variable H i s tr u l y exogeneous. Hence, Var (y-y) = (R + H 2Z 2 2) + Z ^ U 2 + Z 2 2) 9.11 This i s the same expression as that given by Feldstein's formula for a scalar regression model. 2. If H i s known, then Z 1 2 = Z,, = 0. This assumption i s used in the Kalman model for calcul a t i n g the variance of the 1-step innovation. Equation 9.10 then reduces to Var (y-y) = R + H 2Z 2 2 9. 1 2 Z 2 2 i s the mean square error for the state variable x. It is synonymus with "P" in the Kalman F i l t e r notation. 9.5 I l l u s t r a t i o n of the variance formula for the AR(1) The AR(1) model i s used to i l l u s t r a t e the use of the general formula (eqn. 9.10) for the variance of the forecast error. Application of the Kalman F i l t e r to t h i s model was studied in Chapter 6. The state-space formulation considered here corresponds to l e t t i n g the AR c o e f f i c i e n t be the state variable. In practice, flow predictions several time steps ahead are often desired. Because H and x are both unknown and their estimates cannot be assumed to be independent, 121 eqn. 9.10 should be used. The flow at Hope t y p i c a l l y ranges from 1000 m3/s to 7000 m3/s throughout the year. In practice, the AR c o e f f i c i e n t i s of the order of 1.0. Characteristic values for each term of eqn. 9.10 i s examined below: 1 . R Studies from Chapter 5 indicate that the standard deviation of the measurement error is proportional to the flow i t s e l f . In practice, these errors are usually less than 10% (Water Survey of Canada, personal communication). Hence, R i s taken to be .Oly 2. 2 • 2 2 2 This i s given by the Kalman F i l t e r as the MSE of the AR c o e f f i c i e n t , a. As noted e a r l i e r , a i s approximately 1. Estimates of a are obtained at each time step from the Kalman F i l t e r . A conservative estimate of i t s MSE as a percentage of the true value i s about 10%. This results in L2 = .01. 3 . 2 , , This i s the MSE for H. For time t+2, t h i s corresponds to y t + 1 . The MSE(y t + 1) i s given by the Kalman F i l t e r as HPH' + R. Using the above values for P and R, 2,, = .02y 2. 4. 2 1 2 Since fl and x are p o s i t i v e l y correlated, an estimate for Z 1 2 can be obtained by considering the c o r r e l a t i o n between fl and x. The c o r r e l a t i o n i s of the order of .1 122 to 1, hence I 1 2 i s approximately .OOly to . O l y . Each term of eqn. 9.10 i s expressed in terms of the flow magnitude: ( . O l y 2 + . O l y 2 ) + .02y2(1+.01) + .01y (2y + . O l y ) 9.13 Disregarding terms with lower orders of magnitude yie l d s ( . O l y 2 +.0ly 2) + .02y2 + .02y2 9.14 Expression 9.14 indicates that in practice, a l l terms in the general expression are of the same order of magnitude. Therefore, i t i s not j u s t i f i e d to neglect any part of the general formula of eqn. 9.10. Use of the Kalman algorithm to calculate Var (y-y) i s equivalent to using the terms in the f i r s t bracket only of eqn. 9.14. Thus, the actual variance is s i g n i f i c a n t l y larger i f the correct expression is used. 9.6 Conclusions A scalar observation model i s considered in t h i s chapter, y = Hx + e. A general formula for the variance of the forecast error i s developed when H and x are both unknown and t h e i r estimates are correlated with each other. The formula obtained i s : (R + H 2Z 2 2) + I n (x 2 + 2 2 2 ) + Z12 <2Hx + I 1 2 ) 123 The expression consists of three parts: The f i r s t , accounts for the measurement and state estimation errors (R + £ 2 2 ) * The second part r e f l e c t s the uncertainty in the estimate, fl ( I , , ) ; and the f i n a l part accounts for the c o r r e l a t i o n between fl and x ( £ 1 2 ) . Two special cases of the general formula are: 1. F e l d s t e i n 1 s formula (R + Z 2 2H 2) + I,! (x 2 + Z 2 2) This i s used i f fl and x are uncorrelated which i s appropriate i f H i s t r u l y an exogeneous vari a b l e . 2. Kalman F i l t e r formula (R + Z 2 2 H 2 ) This i s appropriate for the 1-step forecast because H i s a known quantity. The importance of t h i s formula i s i l l u s t r a t e d with an application of the AR(1) model to streamflow prediction at Hope. The variance equation given by the Kalman F i l t e r is inappropriate in t h i s case for two reasons: 1. Flow predictions 2-steps ahead or more require knowledge of future flows (unknown H). 2. The estimates fl and x are correlated with each other, as they are both based on past values of y. The p r a c t i c a l i l l u s t r a t i o n indicates that a l l three parts of the general variance expression are of the same order of magnitude, and hence cannot be neglected. In any forecasting s i t u a t i o n , decisions are often made based on the r e l i a b i l i t y 124 of the future predictions. As the r e l i a b i l i t y i s r e f l e c t e d in the standard error of forecast, the correct expression developed in t h i s chapter should be used. 10. SUMMARY AND CONCLUSIONS ARMAX models are often used to describe stochastic processes such as streamflow phenomena in hydrology. These are time series models with exogeneous inputs. Application of the Kalman F i l t e r to these types of models for flood forecasting i s considered in th i s t h e s i s . The Kalman F i l t e r i s a recursive estimation procedure which gives a li n e a r , minimum variance estimator for the state vector at time t. The estimate i s updated at each time step by making use of incoming noise-corrupted observations. The Kalman algorithm, based on a state-space model, accounts for uncertainties both in the states and in the measurements. The performance of thi s estimation technique depends on sa t i s f y i n g the assumptions of the Kalman state-space model. Use of the Kalman F i l t e r in aerospace applications i s very successful because the physical equations and system dynamics are well known. However, such is not the case in streamflow applications as there are many uncertainties associated with th i s stochastic phenomenon. Choice of the "best" ARMAX model with the "proper" noise s t a t i s t i c s i s often an impossible goal to achieve in engineering practice. Nevertheless, the Kalman F i l t e r i s used even when the assumptions of the Kalman model are not completely s a t i s f i e d . The following section summarizes the main contributions of the thesis to the f i e l d of Kalman F i l t e r i n g in streamflow 125 126 forecasting. Subsequent sections give the conclusions for each study. The aim of thi s thesis i s to further the understanding of the Kalman F i l t e r as applied to hydrologic systems. The p r a c t i c a l i t y and the performance of thi s estimation technique i s examined in the context of ARMAX flow models. It should be pointed out that the objective i s not to ide n t i f y the best hydrologic model for a particular phenomenon. No doubt, more physically based models have been proposed and used with much success, for basins with extensive data base necessary as inputs. Recognition of stochastic elements in hydrologic processes has led to the use of ARMAX models in streamflow modelling. It i s for these situations that the application of the Kalman F i l t e r i s considered. In thi s thesis, i t i s assumed that whichever ARMAX model i s chosen, i t i s an adequate description of the underlying process. 10.1 Summary of Thesis Contributions General p r a c t i c a l problems which aris e in the application of the Kalman F i l t e r are described in the l i t e r a t u r e review of Chapter 2. Several problems which frequently occur in hydrologic applications are investigated in t h i s thesis. There are three main contributions of this research, each summarized in the following paragraphs. Advances in the understanding of the Kalman F i l t e r i s made in terms of: 1. f i l t e r ' s s e n s i t i v i t y to input s p e c i f i c a t i o n , and 127 2. p r a c t i c a l i t y of the f i l t e r for d i f f e r e n t state-space formulations of ARMAX models. Correct expression for the mean square error (MSE) of forecast i s derived for the AR(1) model. It i s shown that under appropriate assumptions, the expression reduces to the Kalman equation for the variance of the innovations. This t h i r d contribution is important when autoregressive models are used to forecast flows several time steps ahead. Under the Kalman model, th i s process v i o l a t e s the assumption of known system matrix H, and use of the Kalman equation underestimates the variance of the forecast errors. The problem of input s p e c i f i c a t i o n , for quantities often unknown in practice i s f i r s t investigated. Results indicate that of the four quantities (x 0, P 0, Q, R), only the combined s p e c i f i c a t i o n QR has p r a c t i c a l e f f ects on the observation forecasts. The s e n s i t i v i t y study gives an insight as to how the noise covariances can be spec i f i e d in order to achieve reasonable forecasting performances for the f i l t e r . In addition, the worst s p e c i f i c a t i o n combination i s also noted. The p r a c t i c a l i t y of the Kalman F i l t e r i s i l l u s t r a t e d through three special cases of the ARMAX model. The three models are used to describe flow phenomenon of the Fraser River at Hope; t y p i c a l of basins whose response c h a r a c t e r i s t i c s are constant, or change slowly through time. The generality of the state-space approach allows f l e x i b i l i t y in model formulation. For forecasting purposes, 128 the studies indicate that the most useful formulation i s to write the ARMAX model as the measurement equation with the ARMAX c o e f f i c i e n t s as the state variables. Hence, the Kalman F i l t e r i s used to give recursive estimates of the co e f f i c i e n t s such that flow forecasts are always made with the latest state estimates. Two main advantages of thi s formulation over other possible ones are less data requirements, and robustness of flow forecasts to poor i n i t i a l knowledge of ARMAX c o e f f i c i e n t s . The Kalman F i l t e r i s a 1-step p r e d i c t o r - f i l t e r estimation technique. However, forecasts for several time-steps ahead are required in practice and the f i l t e r i s often used for making these k-step forecasts. In situations where the system matrix H i s unknown, the variance of the forecast error should not be calculated from the Kalman algorithm. A correct expression for this variance i s developed for the univariate AR(1) model. This expression has important consequences in practice because management decisions are often based on the r e l i a b i l i t y of flow predictions which i s indicated by their mean square errors. It i s also shown that in hydrologic applications, a l l terms in the derived expression are of comparable magnitudes. Conclusions of each investigation are given in the following sections. 1 29 10.2 S e n s i t i v i t y Analysis The Kalman algorithm requires input s p e c i f i c a t i o n for the noise covariances and i n i t i a l conditions of the state vector. As these are usually unknown to the hydrologist, the performance of the f i l t e r with respect to misspecification of these inputs (Q, R, x 0, P 0) i s examined. This problem i s studied by formulating an ARMAX model in state-space notation with the model c o e f f i c i e n t s as the state vector. Streamflow data are generated with chosen noise covariances, Q* and R*. Performance of the f i l t e r , based on the observation forecast error, i s examined with respect to input s p e c i f i c a t i o n s . It is found that: 1. For the state-space formulation used, i n i t i a l s p e c i f i c a t i o n of the model c o e f f i c i e n t s are not important. Poor choices of x 0 and P 0 have l i t t l e influence on the flow predictions. 2. Of the four input factors (Q, R, x 0, P 0 ) ; only the combined s p e c i f i c a t i o n of the noise covariances have an important effect on the flow forecasts. 3. If Q and R are unknown, i t i s best to specify them both larger than their expected values. Forecasts obtained are comparable to those for the optimal case of known Q and R. Other combinations of specifying Q and R result in worse forecasting performance than the above. 4. For t h i s state-space formulation, i t i s important to estimate R c o r r e c t l y . The f i l t e r i s i n d i f f e r e n t to 1 30 misspecification of Q i f R i s given c o r r e c t l y . 5. However, the reverse i s not true. Even i f the true Q i s given, the f i l t e r performs worse i f R i s specified too small. Thus, R needs to be estimated also. 10.3 Maximum Likelihood Estimation of the noise variance For each of the three models considered, the method of maximum l i k e l i h o o d i s used to estimate the noise variance. This method i s chosen because i t gives consistent and asymptotic e f f i c i e n t estimates. Evaluation of the log lik e l i h o o d function i s f a c i l i t a t e d by using the Kalman F i l t e r . This i s another application of the Kalman F i l t e r , other than that of forecasting. 10.4 The Autoregressive Model The 1-day ahead forecast for the AR(1) model i s examined. Three schemes are used to forecast the flow at Hope, each corresponding to a possible formulation of the AR(1) model into state-space format. Results indicate that: 1. The best forecasting scheme i s to use the Kalman F i l t e r to obtain updated estimates of the AR c o e f f i c i e n t . Thus, a is the state variable in the state-space framework. Flow predictions at each time step are made with the late s t state estimate. This formulation results in the best forecasting performance. In addition, the flow forecasts are insensitive to i n i t i a l s p e c i f i c a t i o n of a. 131 2. Formulating the AR process as the state equation of the state-space model i s equivalent to not using the Kalman F i l t e r at a l l . An estimate of the AR c o e f f i c i e n t i s required prior to forecasting. Comparable forecasting performance to the best scheme can be obtained for a p a r t i c u l a r value of a. However, the forecasting performance i s sensitive to the s p e c i f i c a t i o n of the AR coef f i c i e n t . 3. The worst performance i s obtained for the scheme which " s p l i t s up" the AR(1) process in the state-space formulation. Flow i s modelled as the state variable while the error term i s modelled as the observation noise. Not only i s the forecasting performance the worst, i t is p a r t i c u l a r l y sensitive to the choice of a determined prior to forecasting. As part of the Kalman F i l t e r algorithm, the 1-step forecast for the observation, and i t s mean square error (MSE) are given. The expressions are y=Hx, and HPH' + R respectively. However, in daily water management, a longer forecasting horizon is usually required. Flow prediction more than 1 day in advance introduces the problem of unknown future flows for autoregressive models. For instance, the 2-step forecast requires knowledge of the flow 1-day ahead. The Kalman algorithm assumes that tomorrow's flow i s known. In practice, a prediction can s t i l l be obtained using the 132 1-step f o r e c a s t from the Kalman F i l t e r as an e s t i m a t e . The MSE f o r the 2-step f o r e c a s t however, cannot be c a l c u l a t e d v i a t h e Kalman a l g o r i t h m . The MSE i s i m p o r t a n t as i t r e f l e c t s t h e h y d r o l o g i s t ' s c o n f i d e n c e i n h i s f o r e c a s t s . Because c o s t l y d e c i s i o n s may depend on the r e l i a b i l i t y of the f l o w p r e d i c t i o n s , the p r o p e r e x p r e s s i o n f o r the v a r i a n c e of t h e f o r e c a s t e r r o r s h o u l d be used. The f i r s t approach t o t h i s problem i s t o c o n s i d e r a d i f f e r e n t ARMAX model which does not r e q u i r e f u t u r e f l o w s i n o r d e r t o p r e d i c t the f l o w 2-days i n advance. T h i s i s d i s c u s s e d i n the next s e c t i o n . The second approach i s t o d e r i v e the c o r r e c t e x p r e s s i o n f o r the mean square e r r o r of f o r e c a s t . T h i s i s summarized i n s e c t i o n 10.7. 10.5 T r a n s f e r F u n c t i o n Model A t r a n s f e r f u n c t i o n model u s i n g upstream f l o w i n p u t s l a g g e d 2 days behind t h a t of Hope i s used. Two f o r e c a s t i n g schemes a r e compared, one w i t h o u t and one w i t h the Kalman F i l t e r . The l a t t e r scheme f o r m u l a t e s the i n p u t - o u t p u t model w i t h t h e c o e f f i c i e n t s as the s t a t e v a r i a b l e s . C o n c l u s i o n s a r e : 1. F o r c o n s t a n t c o e f f i c i e n t s , use of the Kalman F i l t e r does not improve the f l o w f o r e c a s t s . 2. However, f o r e c a s t s o b t a i n e d w i t h o u t the f i l t e r a r e s e n s i t i v e t o the i n i t i a l s p e c i f i c a t i o n of t h e s e c o e f f i c i e n t s . W i t h o u t the f i l t e r , t h e s e have t o be 133 estimated from past data, thus making the choice of dataset important. 3. Although the f i l t e r does not improve the forecasting performance of thi s model, flow predictions are robust to poor s p e c i f i c a t i o n s of the model c o e f f i c i e n t s . Large amounts of data are not necessary for estimating these c o e f f i c i e n t s prior to forecasting. 10.6 Combined ARMAX Model Using a d i f f e r e n t model to forecast k steps ahead for various k's involves too many models as k gets large. Therefore a combined model i s considered in order to achieve a parsimonious representation of the process. This combined ARMAX model i s a combination of the AR(1) and regression models above. The Kalman F i l t e r i s applied to thi s model to give estimates of the ARMAX c o e f f i c i e n t s recursively. Both 1 and 2-step forecasts are obtained and are compared to those of the AR(1) and the transfer function models respectively. Calculation of the 2-step forecast presents the same problem as the AR(1) model, in that tomorrow's flow.is unknown. Nevertheless, the 1-step forecast given by the Kalman F i l t e r i s used as an estimate. The study shows that: 1. The 1-step forecasting performance i s comparable to that of the AR(1) model. 2. The 2-step forecasting performance i s better than that of the transfer function model. Thus, in terms of indentifying the s t a t i s t i c a l model for 1 34 p r e d i c t i n g s t r e a m f l o w s up t o 2 days i n advance, t h i s combined ARMAX model i s adequate. The i m p o r t a n t r e s u l t of t h i s study i s t h a t r e a s o n a b l e f o r e c a s t i n g performance i s o b t a i n e d f o r the 2-step f l o w p r e d i c t i o n s ; even though an e s t i m a t e i s used f o r tomorrow's f l o w i n the system m a t r i x H. V i o l a t i o n of the assumption of known H s t i l l r e s u l t s i n r e a s o n a b l e f l o w p r e d i c t i o n s . 10.7 V a r i a n c e of the f o r e c a s t e r r o r A f o r m u l a f o r the v a r i a n c e of the f o r e c a s t e r r o r when b o t h H and x a r e unknown, i s de v e l o p e d f o r the AR(1) model. I t i s shown t h a t the new e x p r e s s i o n c a l c u l a t e s a v a r i a n c e s i g n i f i c a n t l y l a r g e r than t h a t g i v e n by the Kalman F i l t e r . T h i s e x p r e s s i o n f o r the AR(1) model i s : (R + H 2 Z 2 2 ) + Z,, ( x 2 + Z 2 2 ) + I 1 2 (2Hx + Z, 2) A l l terms i n t h i s e q u a t i o n a re of the same magnitude i n p r a c t i c e . The f i r s t p a r t of t h i s e x p r e s s i o n c o r r e s p o n d s t o the case of known H, and i s e q u i v a l e n t t o the Kalman e q u a t i o n f o r the v a r i a n c e of the i n n o v a t i o n s . Z 2 2 i s the s t a t e e r r o r c o v a r i a n c e m a t r i x . The second p a r t , a c c o u n t s f o r the f a c t t h a t H i s unknown but i t s e s t i m a t e H i s independent of x. Z,, c o n t a i n s the v a r i a n c e s and c o v a r i a n c e s of the elements i n H. The f i n a l term acknowledges p o s s i b l e c o r r e l a t i o n between H and x, as r e f l e c t e d i n Z 1 2 . 135 10.8 State-space Formulation Results of the thesis indicate that application of the Kalman F i l t e r to ARMAX models can give better forecasts when the model c o e f f i c i e n t s are formulated as the state variables. This i s because the f i l t e r updates these c o e f f i c i e n t s recursively in the l i g h t of the observation forecast errors. Although i t presents some conceptual d i f f i c u l t i e s , i t i s necessary to choose between a problem which can be handled in practice and one which i s hard to control . The alternate formulation of allowing the flow variable to be the state, requires estimation of the ARMAX c o e f f i c i e n t s prior to forecasting. This emphasizes the necessity for abundant and good quality data, as forecasting performance i s sensitive to the s p e c i f i c a t i o n of the state t r a n s i t i o n matrix, 10.9 Future Directions In t h i s thesis, hydrologic systems which can be described by constant c o e f f i c i e n t ARMAX models are considered. These are appriopriate when modelling streamflow phenomenon for basins with large drainage areas. Kalman F i l t e r i n g i s applied to these models and the resulting flow predictions are better than or equal to those obtained without the Kalman updating of the model c o e f f i c i e n t s . A useful extension to t h i s would be to investigate the forecasting performance of the f i l t e r for more complex systems whose c h a r a c t e r i s t i c s change s i g n i f i c a n t l y over 136 BIBLIOGRAPHY Akaike, Hirotugu. "Markovian Representation of Stochastic Processes by Canonical Variables." SIAM J . Control, vol. 13, No. 1 (Jan. 1973), 163-73. Benjamin, Jack R., and C. A l l i n C o r n ell. P r o b a b i l i t y , S t a t i s t i c s , and Decision for C i v i l Engineers. New York: McGraw-Hill, 1970. Bolzern, P., M. Ferraro, and G. Fronza. "Adaptive Real-Time Forecast of River Flow-Rates from R a i n f a l l Data." J. of Hydrology, 47 (1980), 251-67. Bowles, David S., and W. J . Grenney. "Estimation of Diffuse Loading of Water Quality Pollutants by Kalman F i l t e r i n g . " In Applications of Kalman F i l t e r to Hydrology, Hydraulics, and Water Resources. Ed. Chao-Lin Chiu. Pittsburg: Univ. of Pittsburg, 1978. Box, G. E. P., and G. M. Jenkins. Time Series Analysis: Forecasting and Control. 2nd ed. San Francisco: Holden Day, 1976. Campbell, W. G., S. Ballou, and C. Slade. Form and Style. 6th ed. Boston: Houghton M i f f l i n Co., 1982. Canada. Environment Canada, Inland Waters Directorate, Water Resources Branch, Water Survey of Canada. Surface Water Data, B r i t i s h Columbia, 1980. Ottawa 1981. Canada. Water Resources Branch, Data Section. Magnetic Tape Records, Ottawa, 1984. Chatfield, C. Analysis of Time Series: An Introduction. 2nd ed. London: Chapman and H a l l , 1980. Chiu, Chao-Lin, and E. 0. Isu. "Application of Kalman F i l t e r to Open Channel Flow P r o f i l e Estimation." In Applications of Kalman F i l t e r to Hydrology, Hydraulics, and Water Resources. Ed. Chao-Lin Chiu. Pittsburg: Univ. of Pittsburg, 1978. Chow, Van Te. "Evolution of Stochastic Hydrology." In Applications of Kalman F i l t e r to Hydrology, Hydraulics, and Water Resources. Ed. Chao-Lin Chiu. Pittsburg: Univ. of Pittsburg, 1978. Cleary, James P., and Hans Levenbach. The Professional Forecaster. Belmont: Wadsworth, Inc., 1982. Constable, T. W., and E. A. McBean. "Kalman F i l t e r Modelling of the Speed River Quality." ASCE-J. of the 137 138 Environmental Engineering D i v i s i o n , v o l . 105, No. EE5 (Oct. 1979), 961-78. Coulthard, John, ed. UBC Tape. UBC Computing Centre, March 1981. Cox, D. R., and E. J. S n e l l . Applied S t a t i s t i c s . London: Chapman and H a l l , 1982. de Jong, Piet . UBC IMSL. UBC Computing Centre, Sept. 1980. , and Phelim P. Boyle. "Monitoring Mortality." J. of Econometrics, 23 (1983), 1-16. , and B. Zehnwirth. "Claims Reserving, State-Space Models and the Kalman F i l t e r . " The Journal of the Institute of Actuaries, v o l . 110, part 1, No. 444 (June 1983), 157-81. Dixon, W. J . , chief ed. BMDP S t a t i s t i c a l Software. Berkeley: Univ. of C a l i f o r n i a Press, 1983. Duong, Nguyen, Rah-Ming L i , and Y. H. Chen. "Adaptive Control of Lock and Dam Gate Openings Using a Kalman F i l t e r for Real-Time I d e n t i f i c a t i o n of Upstream-Downstream Stage Relationship." In Applications of Kalman F i l t e r to Hydrology, Hydraulics, and Water Resources. Ed. Chao-Lin Chiu. Pittsburg: Univ. of Pittsburg, 1978. Feldstein, Martin S. "The Error of Forecast in Econometric Models when the Forecast Period Exogeneous Variables are Stochastic." Econometrica, v o l . 39, No. 1 (Jan. 1971), 55-60. Fox, Daniel J., and Kenneth E. Guire. MIDAS. 3rd ed. Michigan: S t a t i s t i c a l Research Lab. Uhiv. of Michigan, Sept. 1976. Gelb, Arthur, ed. Applied Optimal Estimation. By The Analytic Sciences Corporation. 6th ed. Cambridge: MIT Press, 1980. Greig, Malcolm, ed. UBC ANOVAR. UBC Computing Centre, Oct. 1978. Grigoriu, Mircea. "Prediction and Simulation with Autogressive Models in Hydrology." In Applications of Kalman F i l t e r to Hydrology f Hydraulics, and Water Resources. Ed. Chao-Lin Chiu. Pittsburg: Univ. of Pittsburg, 1978. Halfon, E. "Assessment of Uncertainty in Long Term Prediction made through Mathematical Models." In 139 Cleveland: IFAC, Water and Related Land Resource Systems, 1980. Harrison, P. J., and C. F. Stevens. "Bayesian Forecasting." J . Royal S t a t i s t i c a l Society, B, 38 (1976), 205-47. Harvey, A. C. Time Series Models. Oxford: P. Allan, 1981. , and P. H. J. Todd. "Forecasting Economic Time Series with Structural and Box-Jenkins Models: A Case Study." J. of Business and Economic S t a t i s t i c s , v o l . 1, No. 4 (Oct. 1983), 299-315. Hicks, Charles. Fundamentals Concepts in the Design of Experiments. 2nd ed. New York: Holt, Rinehart and Winston, Inc., 1973. Integrated Software Systems Corporation. T e l l a g r a f . version 4. San Diego,CA: ISSCO, 1981. Jazwinski, A. H. "Adaptive F i l t e r i n g . " Automatica, v o l . 5 (1969), 475-85. . Stochastic Processes and F i l t e r i n g Theory. New York: Academic Press, 1970. Johnson, R. A., and Dean W. Wichern. Applied Multivariate S t a t i s t i c a l Analysis. Englewood C l i f f s , N.J.: Prentice-Hall, 1982. Kahl, Douglas R., and Johannes Ledolter. "A Recursive Kalman F i l t e r Forecasting Approach." Management Science, v o l . 29, No. 11 (Nov. 1983), 1325-333. Kalman, R. E. "A New Approach to Linear F i l t e r i n g and Prediction Problems." Transactions of the ASME - J . of Basic Engineering, (March 1960), 35-45. "A Retrospective af t e r Twenty Years: From the Pure to the Applied." In Applications of Kalman F i l t e r to Hydrology, Hydraulics, and Water Resources. Ed. Chao-Lin Chiu. Pittsburg: Univ. of Pittsburg, 1978. Kendall, M. Time Series. London: G r i f f i n , 1973. Kitanidas, Peter, and Rafael L. Bras. "A Study of C o l l i n e a r i t y and Parameter S t a b i l i t y in Rainfall-Runoff Models: Ridge Regression and Kalman F i l t e r i n g . " In Applications of Kalman F i l t e r to Hydrology, Hydraulics, and Water Resources. Ed. Chao-Lin Chiu. Pittsburg: Univ. of Pittsburg, 1978. • . "Error Indentification in Conceptual Hydrologic Models." In Applications of 140 Kalman F i l t e r to Hydrology, Hydraulics, and Water Resources. Ed. Chao-Lin Chiu. Pittsburg: Univ. of Pittsburg, 1978. Lettenmaier, D. P., and S. J . Burges. "Use of State Estimation Techniques in Water Resources System Modelling." Water Resources B u l l . , 12 (1), (1976), 88-99. L i , Rah-Ming, Nguyen Duong, and Daryl B. Simons. "Application of Kalman F i l t e r for Prediction of State Discharge Relationship in Rivers." In Applications of Kalman F i l t e r to Hydrology, Hydraulics, and Water Resources. Ed. Chao-Lin Chiu. Pittsburg: Univ. of Pittsburg, 1978. Linsley, Ray K. J r . , Max A. Kohler, and Joseph L. H. Paulhus. Hydrology for Engineers. 3rd ed. New York: McGraw-Hill, 1982. Logan, L. A., W. N. Lennox, and T. E. Unny. "A Time-Varying Non-Linear Hydrologic Response Model for Flood Hydrography Estimation in a Noisy Environment." In Applications of Kalman F i l t e r to Hydrology, Hydraulics, and Water Resources. Ed. Chao-Lin Chiu. Pittsburg: Univ. of Pittsburg, 1978. Lyon-Lamb, V i c t o r i a . UBC Tapeasy. UBC Computing Centre, March 1982. McLaughlin, Dennis. "Potential Applications of Kalman F i l t e r i n g Concepts to Groundwater Basin Management." In Applications of Kalman F i l t e r to Hydrology, Hydraulics, and Water Resources. Ed. Chao-Lin Chiu. Pittsburg: Univ. of Pittsburg, 1978. Magnus, Jan R., and H. Neudecker. "The Commutation Matrix: some Properties and Applications." The Annals of S t a t i s t i c s , v o l . 7, No. 2 (1979), 381-94. Mair, Susan. UBC Surface. UBC Computing Centre, Nov. 1982. . UBC Plot. UBC Computing Centre, Nov. 1982. Marsalek, J . "Instrumentation for f i e l d studies of urban runoff." Canada - Ontario Agreement Research Program, Proj. 73-3-R, Res. Rep. No. 42, 82pp. Mehra, R. K. "On the I d e n t i f i c a t i o n of Variances and Adaptive Kalman F i l t e r i n g . " IEEE Transactions on Automatic Control, v o l . AC-15, No. 2 ( A p r i l 1970), 175-84. . "Approaches to Adaptive F i l t e r i n g . " IEEE V 141 Transactions on Automatic Control, (Oct. 1972), 693-98. . "Practical Aspects of Designing Kalman F i l t e r s . " In Applications of Kalman F i l t e r to Hydrology, Hydraulics, and Water Resources. Ed. Chao-Lin Chiu. Pittsburg: Univ. of Pittsburg, 1978. . "Kalman F i l t e r s and their Applications to Forecasting." In TIMS Studies in the Management Sciences. Amsterdam: North-Holland, 1979. . "A Survey of Time Series Modelling and Forecasting Methodology." In Real-Time Forecasting / Control of Water Resource Systems. Ed. E r i c Wood. Oxford: Pergamon Press, 1980. Meinhold, Richard, and Nozer D. Singpurwalla. "Understanding the Kalman F i l t e r . " The American S t a t i s t i c i a n , v o l . 37, No. 2 (May 1983), 123-27. Mendenhall, William. Introduction to Linear Models and the Design and Analysis of Experiments. Belmont: Duxbury Press, 1968. , and Richard L. Scheaffer. Mathematical S t a t i s t i c s with Applications. North Scituate: Duxbury Press, 1973. Moore, John B. WATFIV: Fortran Programming with the WATFIV Compiler. Reston, V i r g i n i a : Reston Publishing Co. Inc., 1 975. Moore, Stephen. "Application of Kalman F i l t e r to Water Quality Studies." In Applications of Kalman F i l t e r to Hydrology, Hydraulics, and Water Resources. Ed. Chao-Lin Chiu. Pittsburg: Univ. of Pittsburg, 1978. Moore, R. J . , and D. A. Jones. "Coupled Bayesian - Kalman F i l t e r Estimtation of Parameters and States of Dynamic Water Quality Model." In Applications of Kalman F i l t e r to Hydrology, Hydraulics, and Water Resources. Ed. Chao-Lin Chiu. Pittsburg: Univ. of Pittsburg, 1978. Nakamura, Masatoshi. "Relationship between Steady State Kalman F i l t e r Gain and Noise Variances." Int. J. Systems Science, v o l . 13, No. 10, (1982), 1153-163. Natale, L., and E. Todini. "A Stable Estimator for Linear Models." Water Resources Research, v o l . 12, No. 4 (Aug. 1976), 667-76. Nico l , Tom, ed. UBC MATRIX. UBC Computing Centre, Mar. 1982. Nightingale, Jon. An Introduction to TEXTFORM. UBC Computing 142 Centre, A p r i l 1984. . Textform Layouts. UBC Computing Centre, A p r i l 1984. # TEXTFORM. UBC Computing Centre, A p r i l 1 984. Patry, G i l l e s G., and Miguel A. Marino. "Time-Variant Nonlinear Functional Runoff Model for Real-Time Forecasting." J. of Hydrology, 66 (1983), 227-44. . Letter to author. Dec. 1983. O'Connell, P. E., and R. T. Clarke. "Adaptive Hydrological Forecasting - A Review." Hydrological Sciences B u l l e t i n , 25, 2 (June 198iT, 179-205. 0'Donovan, Thomas M. Short Term Forecasting. New York: John Wiley & Sons, 1983. Quick, M.C, and A. Pipes. "A combined snowmelt and r a i n f a l l runoff model." Canadian Journal of C i v i l Engineering, Vol.3, No.3, 1976, 449-460. Rao, A. Ramachandra, and R. L. Kashyap. "Causality in Hydrologic Systems." In Cleveland: IFAC, Water and Related Land Resource Systems, 1980. Rodriguez-Iturbe, I., Juan Valdes, and J. M. Velasquez. "Applications of Kalman F i l t e r in Rainfall-Runoff Studies." In Applications of Kalman F i l t e r to Hydrology, Hydraulics, and Water Resources. Ed. Chao-Lin Chiu. Pittsburg: Univ. of Pittsburg, 1978. Sage, A. P., and G. W. Husa. "Adaptive F i l t e r i n g with Unknown Prior S t a t i s t i c s . " Proc. Joint Automat. Control, 1969, 760-69. Schweppe, "Evaluation of Likelihood Functions for Gaussian Signals." IEEE Transactions on Information Theory, 1965, 61-70. Sit t n e r , Walter T. "WMO Project on Intercomparison of Conceptual Models used in Hydrological Forecasting." Hydrological Sciences B u l l e t i n , XXI, 1 (March 1976), 203-222. Stevens, John, Tony Buckland, and Bruce J o l i f f e , eds. UBC Fortran. UBC Computing Centre, Nov. 1981. Szollosi-Nagy, A. "An Adaptive I d e n t i f i c a t i o n and Prediction Algorithm for the Real-Time Forecasting of Hydrological Time Series." Hydrological Sciences B u l l . , XXI, No. 1 143 (March 1976), 163-76. , Ezio Todini, and E r i c Wood. "A State-Space Model for Real-Time Forecasting of Hydrological Time Series." J . of Hydrological Sciences, v o l . 4, No. 1 (1977), 61-76. Todini, E. "Mutually Interactive State-Parameter (MISP) Estimation." In Applications of Kalman F i l t e r to Hydrology, Hydraulics, and Water Resources. Ed. Chao-Lin Chiu. Pittsburg: Univ. of Pittsburg, 1978. , and D. B o u i l l o t . "A Rainfall-Runoff Kalman F i l t e r Model." In System Simulation in Water Resources. Ed. G. C. Vansteenkiste. Amsterdam: North-Holland, 1976. , and J. R. Wallis. "A Real Time Rainfall-Runoff Model for an On-Line Flood Warning System." In Applications of Kalman F i l t e r to Hydrology, Hydraulics, and Water Resources. Ed. Chao-Lin Chiu. Pittsburg: Univ. of Pittsburg, 1978. Valdes, Juan B., J. M. Velasquez, and I. Rodriguez-Iturbe. "Bayesian Discrimination of Hydrologic Forecasting Models Based on the Kalman F i l t e r . " In Applications of Kalman F i l t e r to Hydrology, Hydraulics, and Water Resources. Ed. Chao-Lin Chiu. Pittsburg: Univ. of Pittsburg, 1978. Walpole, Ronald E., and Raymond H. Myers. Probability and S t a t i s t i c s for Engineers and S c i e n t i s t s . 2nd ed. New York: MacMillan, 1978. Wood, E r i c , introd. Real-Time Forecasting / Control of Water Resource Systems. IIASA Proceedings Series. Oxford: Pergamon Press, 1980. Wood, E r i c . "An Application of Kalman F i l t e r i n g to River Flow Forecasting." In Applications of Kalman F i l t e r to Hydrology, Hydraulics, and Water Resources. Ed. Chao-Lin Chiu. Pittsburg: Univ. of Pittsburg, 1978. , and R. K. Mehra. "Model I d e n t i f i c a t i o n and Sampling Rate Analysis for Forecasting the Quality of Plant Intake Water." In Cleveland: IFAC, Water and Related Land Resource Systems, 1980. , and A. Szollosi-Nagy. "An Adaptive Algorithm for Analyzing Short-Term Structural and Parameter Changes in Hydrologic Prediction Models." Water Resources Research, v o l . 14, No. 4 (Aug. 1978), 577-81. Young, P., and P. Whitehead. "A Recursive Approach to Time-Series Analysis for Multi-variable Systems." Int. 1 44 J . C o n t r o l , v o l . 25, No. 3 (1977), 457-82. Z a b l o s k y , J . P a u l . UBC Xerox 9700. UBC Computing C e n t r e , Oct. 1983.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Kalman filter and its application to flow forecasting
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Kalman filter and its application to flow forecasting Ngan, Patricia 1985
pdf
Page Metadata
Item Metadata
Title | Kalman filter and its application to flow forecasting |
Creator |
Ngan, Patricia |
Publisher | University of British Columbia |
Date Issued | 1985 |
Description | The Kalman Filter has been applied to many fields of hydrology, particularly in the area of flood forecasting. This recursive estimation technique is based on a state-space approach which combines model description of a process with data information, and accounts for uncertainties in a hydrologic system. This thesis deals with applications of the Kalman Filter to ARMAX models in the context of streamflow prediction. Implementation of the Kalman Filter requires specification of the noise covariances (Q, R) and initial conditions of the state vector (x₀, P₀). Difficulties arise in streamflow applications because these quantities are often not known. Forecasting performance of the Kalman Filter is examined using synthetic flow data, generated with chosen values for the initial state vector and the noise covariances. An ARMAX model is cast into state-space form with the coefficients as the state vector. Sensitivity of the flow forecasts to specification of x₀, P₀, Q, R, (which may be different from the generation values) is examined. The filter's forecasting performance is mainly affected by the combined specification of Q and R. When both noise covariances are unknown, they should be specified relatively large in order to achieve a reasonable forecasting performance. Specififying Q too small and R too large should be avoided as it results in poor flow forecasts. The filter's performance is also examined using actual flow data from a large river, whose behavior changes slowly with time. Three simple ARMAX models are used for this investigation. Although there are different ways of writing the ARMAX model in state-space form, it is found that the best forecasting scheme is to model the ARMAX coefficients as the state vector. Under this formulation, the Kalman Filter is used to give recursive estimates of the coefficients. Hence flow predictions can be revised at each time step with the latest state estimate. This formulation also has the feature that initial values of the ARMAX coefficients need not be known accurately. The noise variances of each of the three models are estimated by the method of maximum likelihood, whereby the likelihood function is evaluated in terms of the innovations. Analyses of flow data for the stations considered in this thesis, indicate that the variance of the measurement error is proportional to the square of the flow. In practice, flow predictions several time steps in advance are often required. For autoregressive processes, this involves unknown elements in the system matrix H of the Kalman model. The Kalman algorithm underestimates the variance of the forecast error if H and x are both unknown. For the AR(1) model, a general expression for the mean square error of the forecast is developed. It is shown that the formula reduces to the Kalman equation for the case where the system matrix is known. The importance of this formula is realized in forecasting situations where management decisions depend on the reliability of flow predictions, reflected by their mean square errors. |
Subject |
Kalman filtering Stream measurements Flow meters |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-06-23 |
Provider | Vancouver : University of British Columbia Library |
Rights | 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. |
DOI | 10.14288/1.0062573 |
URI | http://hdl.handle.net/2429/25948 |
Degree |
Doctor of Philosophy - PhD |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1985_A1 N52.pdf [ 5.43MB ]
- Metadata
- JSON: 831-1.0062573.json
- JSON-LD: 831-1.0062573-ld.json
- RDF/XML (Pretty): 831-1.0062573-rdf.xml
- RDF/JSON: 831-1.0062573-rdf.json
- Turtle: 831-1.0062573-turtle.txt
- N-Triples: 831-1.0062573-rdf-ntriples.txt
- Original Record: 831-1.0062573-source.json
- Full Text
- 831-1.0062573-fulltext.txt
- Citation
- 831-1.0062573.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0062573/manifest