KALMAN FILTER AND ITS APPLICATION TO FLOW FORECASTING by PATRICIA NGAN B. A. Sc., University of British Columbia, 1981 M. S., California Institute of Technology, 19&2 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in FACULTY OF GRADUATE STUDIES (Civil Engineering) We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA April 1985 Patricia Ngan, 1985 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of 6lQvt- &Nl6t^G:eftlM| The University of British Columbia 1956 Main Mall Vancouver, Canada V6T 1Y3 5-6 f^/R-n ABSTRACT 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 (x0, P0). 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 x0, P0, 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. i i 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. iv Table of Contents Chapter Page ABSTRACT ii LIST OF TABLES x LIST OF FIGURES xiACKNOWLEDGEMENTS x i i i 1. INTRODUCTION 1 1.1 Problem definition 1 1.2 General Objectives 2 1.3 General Thesis Outline 6 2. LITERATURE REVIEW 8 2.1 Explanation of the Kalman Filter 8 2.2 General Applications 12 2.3 Hydrologic Applications 3 2.4 Problems in hydrologic applications 14 2.4.1 Unknown dynamics 12.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 3 3.2 Experimental Procedure 25 3.2.1 Generation of Flow Data 25 3.2.2 Kalman Filter Specifications 27 3.2.3 Factorial Experiment 28 3.2.4 Analysis of Variance 32 v 3.3 -Results 34 3.3.1 ANOVA Results 33.3.2 Two-dimensional Graphs 36 3.3.3 Three-dimensional plot 42 3.4 Discussion of Results 44 3.4.1 ANOVA Results3.4.2 Interpretation of Graphs 45 3.5 Summary 44. HYDROLOGIC SYSTEMS 48 4.1 Time-invariant Systems 49 4.2 Three ARMAX Models 51 4.3 Scope of Applications 3 5. ESTIMATION OF MEASUREMENT NOISE VARIANCE 54 5.1 The Estimation Method 55.2 Estimation Procedure 6 5.3 Transformation of flow data to achieve Uniform Variance 57 5.3.1 Relation between the Noise Variance and the Flow 8 5.3.2 Results 59 5.3.3 Conclusions 60 5.3.4 Box-Cox Transforms 65.4 Maximum Likelihood Estimation of the Noise Var iance • 62 5.4.1 Theory 65.4.2 Results 4 5.5 Summary 5 6. APPLICATION: AR(1) 67 vi 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 coefficient 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 97.1.2 Properties of Scheme 2 0 7.2 General Discussion 91 7.3 Experimental Procedure 92 7.4 Results 93 7.4.1 Estimate of the model coefficients 93 7.4.2 Estimate of the noise variance 94 7.4.3 Performance Indicators for Scheme 1 95 vii 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 definition 119.3 Covariances of Products of normal random variables 117 9.4 Variance of the forecast error for the scalar case 118 9.5 Illustration 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 Sensitivity Analysis 129 10.3 Maximum Likelihood Estimation of the noise variance 130 10.4 The Autoregressive Model 1310.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 13BIBLIOGRAPHY 137 ix LIST OF TABLES Table No. Title Page 3.1 Values of x0*, 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 coefficient 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 76.5b Values of PI2 for Scheme 1 7 6.5c Values of PI3 for Scheme 1 78 6.6 Range of values for a 76.7a Values of PI, for Scheme 2 79 6.7b Values of PI2 for Scheme 26.7c Values of PI3 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 PI2 for Scheme 3 2 6.9c Values of PI3 for Scheme 3 87.1 Least Squares Estimates of the model 94 coefficient 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 Sensitivity of the Performance Indicators 96 to the model coefficient 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 specification 106 8.2b Variation of PI2 to noise specification 106 8.2c Variation of PI3 to noise specification 107 8.3a Variation of PI, to noise specification 107 8.3b Variation of PI2 to noise specification 107 8.3c Variation of PI3 to noise specification 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 Fig. No. Title 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 Filter procedure at each time step. 2.2 Sensitivity 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 fications 4.1 A 'data window' approach to slightly 50 time-invariant hydrologic systems. xii ACKNOWLEDGEMENTS I wish to thank Dr. S.O. (Denis) Russell (Department of Civil Engineering), and Dr. Piet de Jong (Faculty of Commerce) for their valuable guidance and advice throughout this research. The knowledge I have gained is deeply appreciated, and the experience has been most rewarding. Sincere thanks go to Dr. M.C. Quick (Department of Civil Engineering) and Dr. W. Caselton (Department of Civil Enginnering) for their comments on streamflow modelling during this work. Appreciation is 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 requisition. Thanks are also extended to Dr. Glen Cooper (Computing Centre at UBC) and Mr. Gordon Finlay (Department of Civil Engineering, Graphics Lab) for their assistance in the computing aspects of this research. I am grateful to my friends for their informative comments throughtout this thesis. Finally I wish to thank my mother for her support and encouragement, without which this 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. Statistical 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 Filter is an example of an adaptive forecasting scheme. This estimation method is 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 Filter have been very successful in communications and aerospace engineering, fields 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 Filter may be greatly affected when these equations and dynamics are unknown. 1.1 Problem definition Although the Kalman Filter is explained in Chapter 2, it is briefly defined here in order to state the general problem addressed in this 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 is 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 is corrupted by white noise, w.tA>(0,Q) with covariance matrix Q. Similarly, the vector of observations is 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 filtered estimate x is given by the Kalman Filter at every time step. In addition, the mean square error for the state estimate is given. A problem which often arises in practice is that of unknown dynamics; either system matrices H or noise covariances Q , R. As O'Connell and Clarke (1981) point out, this difficulty is often overlooked in applications even though the state estimation procedure depends critically on the proper specification of dynamics. 1.2 General Objectives This research work deals with the problem of unknown dynamics. ARMAX models for describing streamflow phenomenon are considered in this thesis. These are time series models with exogeneous inputs. 3 The practical application used in this thesis is motivated by flood forecasting in B.C. The chosen problem is streamflow prediction for the Fraser River at Hope. This station is located near the downstream portion of the Fraser River, (see Fig. 1.1) As a result 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 typical value of 5000 m3/s. Fig. 1.2 gives an outflow hydrograph for 1983. Station near Spences Bridge Station at Texas Creek /Thompson \ \ River / W. Vancouver' .TOLL CAMtOlA ISIA*0' "A* ""VANCOUVER ^ O Richmond n VAIOU 'O 1 in = 4 5 km = 28 mi Fig. 1.1 Geographic location of the study area. OUTFLOW HYDROGRAPH 8000 I Time Fig. 1. 2 Outflow hydrograph at Hope, B. 6 In this thesis, ARMAX models are cast into the state-space format with the coefficients as the state vector. In hydrologic applications, future flow is the quantity of interest. Forecasting performance of the filter is measured in terms of the observation forecast error (y-y). The general objectives of the thesis are: 1. to examine the sensitivity of the Kalman Filter with respect to mis-specifications of noise covariances and initial conditions of the state vector. These inputs are required by the Kalman algorithm, but are often unknown in practice. 2. to investigate the maximum likelihood method for estimating the noise covariances in hydrologic models of practical interest. 3. to compare the forecasting performance of ARMAX models depending on whether the Kalman Filter is used or not. 1.3 General Thesis Outline The next chapter reviews the literature on the Kalman Filter. Included in this review are applications and the problems which arise in practice. Chapter 3 investigates the sensitivity of the filter 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 likelihood 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 Filter is used. Often in forecasting, the standard deviation of the prediction is required; but there are many cases where the expression for the variance given by the Kalman Filter is not applicable. Chapter 9 discusses this in detail and a theoretical expression is presented for the variance of the forecast error. The importance of this is illustrated in an example. Finally, chapter 10 summarizes the results of this thesis and gives recommendations for future research. 2. LITERATURE REVIEW There are five sections to this chapter. First, the Kalman Filter is described in detail. General applications are given in section 2 and hydrologic examples in section 3. Practical problems in the implementation of the filter are discussed in section 4. Finally, section 5 outlines in detail the specific topics addressed in subsequent chapters. 2.1 Explanation of the Kalman Filter The following notation is used in this thesis. Vectors are denoted by underlined, small letters; and matrices are denoted by capital letters. The linear state-space model, sometimes known as the Gauss-Markov model consists of two equations: State eqn: xfc = $ Xfi + 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 is given by 3>, the state transition matrix. The states are corrupted by white noise, w distributed with mean 0 and covariance matrix Q. $ and Q can be time-varying, but are assumed known at all 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 distributed with mean 0 and covariance matrix R. Similarly, H and R may be time-varying, but are assumed known at all times. Moreover, wfc and vfc are assumed to be serially and mutually uncorrelated. For the Kalman model consisting of eqn. 2.1 and 2.2, the objective is to estimate the state vector at current time, xfc given all past measurements. Obtaining an estimate for the state at the current time, is known as filtering. In particular, if all past information can be summarized in a prior estimate for the state vector, this type of estimation is known as recursive filtering. The Kalman Filter is a recursive procedure which yields estimates for xfc at every time step. This filtering 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, it is not a Bayesian technique per se, as often quoted in the literature, despite this decision theoretic approach. In fact, xt (the 'posterior' or updated estimate) is derived from a least squares criterion by minimizing the residual of the state estimation vector. The Kalman Filter is a solution which yields a linear, minimum variance estimator for the state vector at every time step. In addition, the error covariance matrix P , associated with x., is given. 10 The Kalman Filter algorithm can be divided into two pa r t s: 1 . Prior to observing y_ -t/t-1 = *t-1 -t-1 2,3 Pt/t-1 = *t-1 Pt-1 *'t-1 + Qt 2'4 2. After observing y_t, it = h/t-i + Kt[*t - £t ] 2-5 Pt = [I - Kfc Ht] 2.6 where Kfc = P^^ H't i^t/t^»\ + *t ]' 1 2.7 The 1-step forecast for the observation is: £t+1 = HA+l/t 2'8 The Kalman Filter accounts for uncertainties by providing an algorithm which sequentially combines model (state equation) and data (measurement equation) information to yield updated estimates of the state vector. These estimates can be projected forward to obtain future predictions of the observations. It is computationally appealing due to its recursive nature, i.e. all previous information is contained in the prior estimate of the state. A schematic representation of this procedure is shown in Fig. 2.1. 11 State equation at time t-1 1 ~ Q Observations J • at time t R i Fig. 2.1 A Schematic representation of the Kalman Filter procedure at each time step. 12 2.2 General Applications Application of the Kalman Filter is widespread, and the literature spans over many disciplines. It is used particularly in forecasting situations. The classic example is from aerospace engineering. The problem is to continually estimate the position, velocity and acceleration of a target in Cartesian co-ordinates, upon receiving noise corrupted measurements from a radar in polar co-ordinates. In finance, Kalman Filtering has been applied to estimate the regression coefficients 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 this are found in the work by Young and Whitehead (1977), Constable and McBean (1979). However, non-linearities in these formulations lead to the use of the extended Kalman Filter which is not the subject of the thesis. The Kalman Filter has also been illustrated in many aspects of the hydrologic process. Bras (1978) has successfully applied this technique to sampling network design. An example of estimating groundwater basin characteristics is given by McLaughlin (1978). However, the dimensionality of groundwater models is large, as many parameters need to be estimated. Computational efficiency and the problem of joint state and parameter estimation are two major difficulties (Wilson et al., 1978). This filtering technique has also been used in the operation of water resource systems. An example is the 13 maintaining of water levels for navigation purposes. Duong et al.(l978) applied the Kalman Filter to obtain real time estimates of system parameters required for controlling a lock and dam gate. 2.3 Hydrologic Applications The most abundant applications of this estimation technique are in streamflow forecasting. In many situations, the underlying process can be described by a statistical model such as a regression or time series model. In these cases, the Kalman Filter is a useful supplement to hydrologic forecasting. The rainfall-runoff process is often modelled by an autoregressive scheme plus some input information. There are two ways of formulating this type of hydrologic model into the state-space notation. Rodriguez-Iturbe (1978) has formulated the problem with the coefficients of an ARMA model as the variables of interest. This is also shown in an example by Wood et al., (1978). The instantaneous unit hydrograph, IUH represented by a convolution integral, is a classic example of such a forecasting scheme. For examples, see Rodriguez-Iturbe (1978), Szollosi-Nagy (1976). Alternatively, the discharges can be taken as the states, but the unknown ARMAX parameters then need to be estimated a priori. An example of this is given in Szollosi-Nagy et al.,(l977). Conceptual response models involving the continuity and storage discharge equations have also been 14 used with the Kalman Filter. However, these models lead to non-linear estimation. Wood (1978) investigated flood routing models via the Kalman Filter for forecasting water levels. 2.4 Problems in hydrologic applications 2.4.1 Unknown dynamics The problem of unknown noise covariances is often encountered in hydrologic modelling as a result of the complexity of the runoff process. Mehra (1979) pointed out the extreme sensitivity of the Kalman Filter 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 in measurement noise variance on steady-state RMS error in state esti mation. or-= measurement noise standard deviation, (from Mehra 1968). J_ attumad/O. actuat q q Effect of errors in process noise variance on steady-state RMS error in state estimation. o°n = process noise standard deviation, (from Mehra 1968). Fig. 2.2 Sensitivity of the state estimates to mis-specifications of noise standard deviations. 15 This was shown in qualitative terms and the idea conforms with one's intuition. However, the filter's performance under simultaneous misspecification of Q and R needs to be examined. In practice, these noise covariances are specified too small or too large with respect to their true values. The manner in which these specifications affect the filter's performance is important knowledge for the forecaster. Moreover, the initial conditions of the state are often unknown. These questions are addressed in the next chapter. Because flow is the variable of interest, the filter'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. $ is the state system matrix. When elements of # are unknown, this leads to non-linear filtering. One way of resolving this is to linearize the problem with respect to the latest state estimate by the Extended Kalman Filter. However, Mehra has indicated that the estimates thus obtained are sensitive to initial conditions. Moreover, the EKF has not been proven to yield 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 non-linearities of the ARMA coefficients in An extension to this 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 Filters: one which yields the state estimates (given the parameter values) and the other which yields parameter estimates (given the state values). A constrained linear estimator proposed by Natale and Todini (1976) was illustrated 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, this necessarily leads to a different expression for the variance of the forecast error. Authors have overlooked this (Harrison & Stevens, 1976) when they suggested using the same Kalman Filter formula to obtain the variance. Whether or not this is a serious error in practice depends on the application. Feldstein (1971) has studied this problem in the context of econometrics. However, because one of his assumptions is often not satisfied in hydrologic applications, this problem needs further investigation. 2.4.4 Estimation of the noise covariances Implementation of the Kalman Filter allows the forecaster to exercise his judgement regarding the accuracy of the underlying model vs. that of the observations. This is achieved by specifying the noise covariances, Q and R at each time step. The Kalman gain matrix K, in effect, 17 controls how the state estimates are updated. Hence the performance of the filter depends critically on these values. These are often unknown in water resource applications, and many authors have addressed this problem. In real-time applications where historical data are scarce, adaptive algorithms are often used (Mehra 1973, IEEE). Jazwinski & Bailie (1969), Sage & Husa (1969) propose covariance matching techniques, where the theoretical covariance of the innovations and the sample covariance are matched by adjusting Q and R appropriately. This heuristic approach is computationally attractive though it is not guaranteed to converge. The special case of known Q but unknown R is reported to have been handled more successfully (Mehra, 1972). Correlation methods only apply to time-invariant systems under stationary conditions which impose a limitation for some water resource systems. The Bayesian method is used to compute the posterior probability of the 6 = {Q, R} given the observed data, i.e. p(0|yfc). Calculation of this requires the likelihood of 6, £(fl|yt) = p(yt| yt_.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 Filters. Posterior probability is calculated for each one of these Kalman Filters. As it can be readily seen, this is not attractive computationally unless N and p are both small. This procedure discriminates the best model out of a given group 18 of models. However, if the original set of 0's do not encompass the true value, then the "best model" is still not optimal in the sense that there is still a better 6 that could be used. This method was used by Moore & Jones (1978), and Valdes et al., (1978). The tremendous complexity of the system renders it computationally unattractive and cumbersome. The previously described methods are somewhat ad hoc in their applications. The Maximum likelihood method for estimating Q and R is based on a mathematical principle and can be universally applied. The main advantage of this method is that it is asymptotically unbiased, consistent and efficient. The idea is to maximize the likelihood of observing a particular combination of Q and R given the observation y over a range of values for the noise covariances. The major drawback is that the likelihood expression involves computation with a non-diagonal matrix. Because the previously described methods are not guaranteed to work, this method is investigated in chapter 5. It should be noted that although the problem of estimating noise variances has been addressed theoretically in the literature, the results have had limited use in practice. The main reasons are: 1. Many of the proposed solutions were developed for some particular application. Hence they tend to be application dependent. 2. The methods which are easily implemented computationally 19 are usually heuristic in nature. They have not been proven to give reliable estimates. Moreover, different authors, who have used them, reported contradictory results (Mehra, 1972; O'Connell & Clarke, 1981; Sage & Husa, 1969; Wood et^ al., 1978). 3. The more complex methods are not easily implemented, though in general, they yield more consistent estimates (Mehra, 1972,1980). in addition, the complex mathematics required do not lend themselves to practical implementations. In fact, they act more like a black-box to the forecaster. There is a definite need to bridge the gap between the theoretical research and the implementation feasibility required in practice. 2.5 Detailed thesis objectives and outline As discussed in section 2.4, the noise covariances and the initial specifications for the state vector are often unknown. These correspond to Q, R, x0, Po of the Kalman state-space model. Two questions are addressed in chapter 3. 1. Which of the above quantities (Q, R, x0, P0) if misspecified, have a practical effect on the filter's performance. 2. How do these quantities affect the forecasting performance. The second half of the thesis is based on ARMAX models. Statistical models are chosen for the applications because 20 they are amenable to state-space formulation of the Kalman Filter. Conceptual models such as SWMM and HEC I are not considered in this thesis. These are physically based models and are more difficult to implement. These are more worthwhile to consider for long-term forecasting for a river basin. On the other hand, ARMAX models require relatively less data. This thesis is focused on the Kalman Filter rather than model identification of hydrologic processes. Thus, statistical models are used. Moreover, the simpler representation of streamflow phenomenon allows a more direct investigation of the practical problems associated with the Kalman Filter. Chapter 4 describes three ARMAX models to be used for the remainder of the thesis. The method of maximum likelihood is used to estimate the noise variance for these models in chapter 5. In particular, a simplified approach to evaluating the log-likelihood expression is examined. This has been proposed by Schweppe and illustrated by Ledolter ett al. (1983). The next three chapters investigate the forecasting performance of the three models depending on whether or not the Kalman Filter is 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. Finally, a general 21 expression for the mean square error of the forecast when H is unknown, is developed in Chapter 9. This situation occurs, for example, when the AR(1) model is used to predict the flow than more two time steps ahead. In this case, the variance of the forecast error is not given by the Kalman equation, HPH' + R. 3. SENSITIVITY ANALYSIS The Kalman Filter has been used in many hydrologic applications. It is a popular estimation technique due to its recursive nature and its ability to handle uncertainties. However, its applications have not always been successful. This is often because the assumption of known input quantities, is not satisfied. These quantities are the initial conditions of the state vector x0 , P0, and the noise covariance matrices Q, R. The Kalman algorithm assumes that x0t Poi Q» R are correctly specified. But when they are unknown in practice, their mis-specification may result in unreliable estimates. The objective of this study is to examine the sensitivity of the filter to mis-specification of the quantities x0, P0, Q, R. The sensitivity study is analyzed as a statistical factorial experiment. Each of the input quantities {x0, P0, Q, R} is treated as a factor. The forecasting performance of the filter is measured by an indicator based on the observation forecast error. This indicator is the response variable in the factorial exper iment. A scalar ARMAX model is used for streamflow modelling in this study. In this application, an autoregressive model of order 1, with temperature as an exogeneous variable is chosen. The model is cast into state-space format with the model coefficients as the state vector. Hence Q represents the noise covariance matrix for the coefficients, and R is 22 23 the variance of the measurement error. Results of this sensitivity study show that initial specifications of the ARMAX coefficients, x0 and P0, do not affect the forecasting performance of the filter to any degree of practical concern. The performance is materially affected by the combined specification 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 relatively large results in a filter performance comparable to the case of known Q and R. Specifying Q small and R large yields the worst forecasting performance. If the measurement noise variance is known, then the forecasting performance is indifferent to under or over-specification of Q. However, if Q is known, it is found that better filter performance is obtained if the measurement noise variance is specified relatively large. 3.1 Experimental Plan The approach to this investigation is divided into four parts. 1. Flow generation Flow sequences are generated by a chosen ARMAX model with known input quantities: x0*, Q*, R*. 2. Kalman Filter Specification The flow sequences are treated as observations and a Kalman Filter is applied to each flow sequence. The 24 state-space format used corresponds to allowing the ARMAX coefficients be the state vector in the Kalman model. Hence, estimates of the coefficients are given at each time step and flow forecasts are made with the latest estimate given by the filter. The noise covariance matrices used in the generation of streamflow data are Q* and R*. The Kalman algorithm uses x0, P0, Q, and R where these may be different from the generation values. Q and R can be specified to be less than, equal to, or greater than their true values, Q* and R*. Hence, the sensitivity of the filter's performance to these specifications is examined. The performance of the filter is measured by an indicator based on the observation forecast error. A Kalman Filter is applied to each combination of the input quantities. One of these filters contain the correct specifications of the inputs. Factorial Experiment The investigation is studied as a statistical factorial experiment with four factors; namely {x0, P0, Q, R}. The levels of these factors correspond to the different specifications. Analysis of Variance Sensitivity of the filter's performance is investigated through analysis of variance (ANOVA) on the performance indicators. The analyses indicate which factors or combinations thereof, have a significant effect on the 25 performance indicator. 3.2 Experimental Procedure 3.2.1 Generation of Flow Data Application of the Kalman Filter in this thesis is motivated by real-time flow forecasting of the Fraser River. A stochastic model for describing the streamflow phenomenon is chosen: qt = at *t-1 + bt Tempt + vt 3.1 The autoregressive term is characteristic of stations with large drainage areas, as the flow is dependent on upstream storage effects. The temperature term represents the snowmelt influence on the runoff. The coefficients of the model, afc and bfc, are chosen to follow a random walk. The ARMAX model of eqn. 3.1 is recast into state-space framework as follows: State eqn: V "at-r + "w,t" bt _bt-1_ _w2t. 3.2 Obs. eqn: [qt_1 Tempt] 3.3 The noise terms are chosen with the following properties: 26 1• wfc ~ (0,Q*) The covariance matrix for the coefficients (states) is constant with respect to time. It is also a diagonal matrix. Thus, independent errors for the model coefficients are assumed. 2. vt<v (0,R*) Similarly, 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 particular combination of x_o*, Q* i R* • The following table gives the values used in the flow generation. Table 3.1 Values of x0*, Q*, R* used in flow generation x0 * = .85 7.00 Arbitrary starting values are used as this does not affect the streamflow generation. Qi* = 000036 0 0 .0016 0 * = yII .0001 0 0 .01 0 * = yIII .0004 0 0 .0225 R:* = 100 Rn* = 400 = 625 RIV* = 900 The chosen values of Q* are representative of practical situations; a standard deviation of .5% and 2% from one time 27 step to the next. Similarly, the range of R* reflects a measurement error with a standard deviation from 3% to 12% of the flow. In practice, the standard deviation of the measurement error is not more than 10% of the flow. The number of combinations of the above input quantities is twelve. Hence, the number of CODES is 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 * yII CODE 5 CODE 6 CODE 7 CODE 8 0 * UIII CODE 9 CODE 10 CODE 11 CODE 12 3.2.2 Kalman Filter Specifications Values of the quantities used in data generation are denoted with a '*'; otherwise the letters represent specifications of these quantities for the Kalman algorithm. Therefore, Q* is used in generation; while Q is a specification value. Different specifications for {x0i Por Qf Rl result in different Kalman Filters. For each 28 code, 36 Kalman Filters 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 Filter is measured in terms of the flow forecast error. The indicator is relative root mean square (RRMS) error defined as: This is preferred over the common RMS error as RRMS does not depend on the units of measurement. 3.2.3 Factorial Experiment The investigation is set up as a factorial experiment with a completely randomized design. It is a fully-crossed four factor experiment. The factors are the specifications of the input quantities: {x0r po> Q> R) • The number of levels for each factor is given below: Table 3.3 Number of Levels for the Factors RRMS 3.4 Q R 2 2 3 3 l 29 The levels correspond to different specifications of the input quantities. The number of ways of combining these levels is: k = II / . = 36. Hence, 36 different Kalman i Filters are used under each CODE. The levels of these factors are described below. Table 3.4 Levels of the Noise Covariances level x0 P0 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 x0 corresponds to an initial specification error of 7%, while the bad estimate corresponds to 50%. These limits are representative of the range of accuracies for estimates of x0 in practice. Level 1 of the noise specification corresponds to under-estimating the standard deviation by a factor of 2. Similarly, level 3 30 corresponds to over-estimating by a factor of 2. Level 2 represents a specification which is equal to the true covariance matrix. In practice, it is expected that estimates of Q and R, are within the bounds provided by levels 1 and 3. The response variable of the factorial experiment is 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 effect on the performance of the Kalman Filter. Each of the 36 Kalman Filters is subjected to n=5 series of observations, the 5 series being different in each case. This conforms to a completely randomized (CR) design. For this non-repeated measures design, the measurement errors are uncorrelated. Thus, an unrepresentative series does not replicate its characteristic from one Kalman Filter to another. Alternatively, the same 5 series of observations could be used for each Kalman Filter in order to minimize the variation between groups. This repeated measures design results in correlated measurement errors among the various Kalman Filters. Hence, a CR design is chosen. The statistical model for a factorial experiment with 4 factors is: Y. i j klm = M + A. + B.+C, + D ilk + 2-way interactions .. + 3-way interactions ... + 4-way interaction + e m 3.5 31 For this sensitivity study, the following variable substitution is made: RRMS = M + X0 + P0+Q + R x_oPo + • • • + x0P0Q + ... + x0P0QR + e The notation is described as follows: 1. RRMS is the response variable of the experiment. 2. ju is the overall mean of the performance indicator, RRMS. 3. x0, P0, Q, R are the factors and they are called the main effects in eqn. 3.6 4. The 2-way interactions include all possible combinations of the factors in groups of two. These are x0P0, x0Q, x0R, P0Q, P0R, QR. 5. Similarly, the 3-way interactions are all combinations of the factors in groups of three. These are x0P0Q, x0P0R, P0QR, x0QR. 6. Finally, there is the 4-way interaction and the error, 3.6 32 3.2.4 Analysis of Variance The experiment is analyzed by ANOVA with all terms (except M) on the right hand side of eqn. 3.6 as the 'sources' in the ANOVA table. One analysis is done for each CODE; hence, twelve ANOVA's are made in total. The computer package 'ANOVAR' is used and a sample output of the ANOVA table is 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 qualitative guide for pointing out those sources which are important in practical situations. The p-values are compared among the sources in order to indicate those which have a significant effect on the response, RRMS. The smaller the p-value, the more likely it is that the corresponding source has a significant effect on the filter's performance. A flow chart of the procedure is given in Fig. 3.1. 33 Generate flows choose x R'" o generate N=kn=180 flow sequences of 100 time steps long i Kalman Filter .KF1,1 KF1,2 •• • KF1,5 Specifications KF36,1 KF36,2 " • KF36,5 Factorial Experiment measure the performance of the Kalman Filter, RRMS ANOVA determine significant sources for RRMS. Repeat for each CODE Fig. 3 . 1 Schematic diagram of Experimental Procedure. 34 3.3 Results 3.3.1 ANOVA Results The ANOVA results indicate that under all CODES, the significant sources are the specification of the main effect Q, and the interaction QR. The significance of QR means that it is the combined specification of the noise covariances which affects the filter's performance. A sample output of the ANOVA table for CODE 7 is 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 all CODES. The manner in which the QR interaction affects RRMS is obtained from the following two and three-dimensional graphs. The mean RRMS values are plotted with respect to the levels of the source, QR. 35 Table 3.5 Computer Output ANOVA Results for CODE 7 Factors: A = *0 B = P 0 C = Q D = R p-values; given 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 that the p-values for sources C=Q and CD=QR are .0117 and .0151 respectively. 36 3.3.2 Two-dimensional Graphs The ordinate in the following graphs is the performance indicator, RRMS. This indicator is such that the smaller the value, the better is the performance. The abscissa denotes the three levels of specification 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 level of R (i.e. holding R constant). Those in series (b) are with respect to R for each level of Q. The features displayed by the graphs are summarized below: 1. From series (b) of Figs. 3.2-3.4, it can be seen that the smallest RRMS value occurs at Q(2)R(2). This represents correct specifications of the noise covariances. 2. The graphs in series (a) show that if R is given correctly, then specifying Q too small Q(1), or too large Q(3), results in almost identical RRMS values. In other words, the filter's performance is indifferent to mis-specification of Q if R is known. 3. A different behavior is observed by examining the graphs in series (b) of Figs. 3.2-3.4 for level Q(2). It is noted that if Q is given correctly, then under-specification of R results in a larger RRMS than over-specification. This phenomenon is particularly noticeable under Q*III shown in Fig. 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 relatively close to the minimum RRMS. These values correspond to the specifications, Q(2)R(3) and Q(3)R(3). It is noted that both these specifications involve R(3). This feature is more pronounced as Q* and R* increase. In general, specifying R too small should be avoided. High values of the performance indicator for all specifications of Q under level R(1) are obtained. A large increase in the RRMS value is noted for Q(1)R(3). For example, see Fig. 3.3 (b). 38 Fig. 3.2 (a) RRMS with respect to Q Fig. 3.2 (b) RRMS with respect-to R 39 CODE 5 CODE 5 0.03*. 0.03-t Legend 0.08 CODE 8 Legend o ot » Q2_ & 03 Fig. 3.3 (a) RRMS with respect to Q Fig. 3.3 (b) RRMS with respect to R 40 Fig. 3.4 (a) RRMS with respect to Q Fig. 3.4 (b) RRMS with respect to R Fig. 3.4 (a) cont'd RRMS with respect to Q Fig. 3.4 (b) cont'd RRMS with respect to R 42 3.3.3 Three-dimensional plot A three-dimensional plot of the interpolated response surface for the RRMS is given. This is plotted with respect to the nine specification levels of both Q and R. Because the behavior of the filter to QR interaction is the same, only one response surface is drawn. This is a qualitative 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 rises in all directions away from the minimum. It is noted that the topography is not too different 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 specification, Q(1)R(3). Q levels RRMS R levels A3 Q(1)R(3) Fig. 3.5 Response Surface RRMS with respect to QR specifications 44 3.4 Discussion of Results 3.4.1 ANOVA Results The factorial experiment considers four factors: x0, P0, Q, R. It is found that: 1. Initial specifications of the ARMAX coefficients x0, P0 do not have a significant effect on the filter's forecasting performance. In practice, these coefficients are often unknown and traditionally, much work has gone into their estimation prior to forecasting. This study shows that even an error of 50% in their initial specification, does not significantly degrade the performance. 2. With the coefficients as the state variables in the Kalman model, it is found that the filter is sensitive to the combined specification of the noise covariances, Q and R. The 1-step flow forecast is based on the latest 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 ratio of the uncertainties in the state estimate to the measurement noise. Because the uncertainty in the state estimate is affected by the noise covariance Q, Kfc can be considered as a ratio of Q to R. Hence, it is the combined specification of the noise covariances which affects the filter's performance. 45 3.4.2 Interpretation of Graphs Because of the interaction between Q and R, it is not sufficient to examine the filter's sensitivity to either Q or R individually. It should be examined with respect to both quantities. 1. The optimum level of performance (represented by the minimum RRMS value) is obtained when Q and R are specified correctly. This is true under all testing conditions. 2. If the variance of the measurement error R is known, then the filter's performance is indifferent to under or overspecification of Q. Its performance level can only be improved if Q is specified correctly. 3. If Q is known to begin with, then specifying R too small should be avoided. 4. For systems with large disturbances in the ARMAX coefficients, under-specification of Q should be avoided. This is particularly true when R is specified large. The worst performance occurs for Q(1)R(3). 5. A general guideline is to specify both covariance matrices large because relatively good forecasting performance can be obtained under Q(3)R(3). 3.5 Summary The Kalman Filter is applied to an ARMAX model for streamflow forecasting. A common problem is that initial conditions of the model coefficients and the noise 46 covariances are unknown. The sensitivity of the filter's forecasting performance to the specification of x 0, P0, Q, R is examined. The ARMAX model is recast into state-space format with the coefficients as the states in the Kalman model. For this sensitivity study, the system matrices $ and H of the Kalman Filter are known. The filter's performance is measured in terms of the observation forecast error, defined by RRMS. The - conditions tested correspond to a standard deviation of .5 to 2% change in the model coefficients, and a standard deviation of 3 to 12% in the measurement noise. It is found that initial specifications of the ARMAX coefficients do not affect the filter's performance. This feature is useful for the practicing engineer; especially when there is insufficient data for obtaining reliable estimates of the coefficients prior to forecasting. The study also shows that the filter's performance is sensitive to the combined specification of the noise covariances. If Q and R are unknown, it is best to specify both covariances larger than their expected values. This specification results in a forecasting performance comparable to the case of known Q and R. If only one of the noise covariances can be estimated, it is better to estimate R accurately. This is due to the filter's insensitivity to under or over specification of Q, if R is given correctly. However, if Q is given correctly, the filter performs worse if R is specified too small. 47 In general, it is best to specify a large variance for the measurement error. This results in good forecasting performance as long as Q is not specified much smaller than its true value. 4. HYDROLOGIC SYSTEMS The Kalman Filter has been applied in the previous sensitivity study to a general ARMAX model, used to generate synthetic streamflow data. The state-space formulation corresponds to letting the ARMAX coefficients be the state variables in the Kalman model. The study shows how the forecasting performance of the filter is affected by the combined noise covariance specification. It also indicates that initial specification of the ARMAX coefficients do not affect the flow forecasting performance. This chapter describes three particular ARMAX models, chosen for Kalman Filter applications in Chapters 6 to 8. Actual daily flow data is used in these studies for forecasting flows 1 and 2 days in advance. These models represent hydrologic systems with certain basin characteristics, described in this chapter. The practicality of the Kalman Filter cannot be explored for all types of hydrologic systems represented by the general ARMAX model, as it is data dependent to some extent. Hence, specialization to certain types of systems is necessary. The practicality of the Kalman Filter and its performance under different state-space formulations of a hydrologic model are investigated using real data. Studies in Chapters 6 to 8 use daily flow data from the Fraser River at Hope in B.C. The Fraser River at Hope is one of the largest rivers in Canada and the flow changes relatively slowly through time. As discussed in section 4.1, 48 49 this allows the assumption of time-invariant systems. Therefore, ARMAX models with constant coefficients will 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 relatively 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 river 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 is on the Kalman Filter and its 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 practical interest in water resource management. 4.1 Time-invariant Systems Flow forecasting using ARMAX models with time-invariant coefficients are considered. Constant coefficient 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 effects. For these systems, rainfall inputs have little influence on the outflow. Rodriguez-Iturbe (1978) and Szollosi-Nagy (1976) concluded that the assumption of time-invariant coefficients can be extended to systems whose response characteristics change slowly over time. A "window" type of approach is suggested for handling the time varying nature, as shown in the diagram. TIME (from Szollosi-Nagy, 1976) z(t) is the observation at time t. Fig. 4'. 1 A 'data window' approach to slightly time-variant hydrologic systems. 51 Application of the Kalman Filter in this thesis is motivated by the following problem. Streamflow predictions for the Fraser River at Hope are required in the control of this water resource system. The maximum daily flow of the year (peak runoff) which occurs in June, is mainly due to snowmelt. Therefore, a reasonable forecasting horizon is April 1 to September 30. Within this time frame, it is expected that the basin characteristics do not change abruptly. The station at Hope is located near the downstream end of the Fraser River. This means that it has a large drainage area and thus a slow basin response. Hence, time-invariant ARMAX models are appropriate for the investigations in this thesis. As a result, only one noise variance is to be estimated, and this is addressed in Chapter 5. 4.2 Three ARMAX Models Three stochastic models commonly used in streamflow modelling are considered. 1. AR(1) qt = aqt_1 + vfc 2. TRANSFER FUNCTION qt = biq*t-2 + b2<3**t-2 + Vt 3. Combined, ARMAX qt = c1qt_1 + c2q*t_2 + c3q**t_2 + vfc where q, q*, q** denote streamflows at three different 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 justifications 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 effects can be represented by autoregressive terms. Predictions beyond this can be made if the unknown flows at lag one are replaced by their estimate. However, this necessarily results in a larger variance for the flow forecast error. This topic is addressed in Chapter 9. Application of the Kalman Filter to this model is investigated in Chapter 6. 2. Upstream inputs can contribute to the rising portion of an outflow hydrograph as pointed out by Wood (1978). A transfer function model can be used where the coefficients b, and b2 represent flow contributions from upstream stations. The 2-day lag in the flows allow the prediction of the 2-day forecast with all the exogeneous variables known at time t. The application of the Kalman Filter to this model is examined in Chapter 7. In fact, flood routing procedures described by the convolution integral are often linearized. The discrete time representation of the integral can be considered as a transfer function model. An example of this 53 formulation can be found in Wood(l978). 3. The third model combines the autoregressive process with upstream inputs. The Kalman Filter is applied to this model in Chapter 8. 4.3 Scope of Applications The state-space formulation whereby the ARMAX coefficients form the elements of the state vector is considered. This format is used throughout the remaining chapters, except where otherwise noted. The problems associated with the specification of the inputs {x0, P0, Q, R} have been reduced to the problem of estimating R only because: 1. Q=0 as time-invariant coefficients are considered. 2. Initial specification of the ARMAX coefficients, x0, P0 do not affect the filter's forecasting performance. 5. ESTIMATION OF MEASUREMENT NOISE VARIANCE Results of Chapter 3 show that the filter's forecasting performance is affected by the specification of the noise covariance matrices. In fact, variance estimation in state-space models is a subject of much current interest. A review of this was given in Chapter 2, section 2.4.4. In this chapter, a method for estimating the noise variance is investigated. The three models discussed in Chapter 4 are considered. It is found that the variance of the measurement error cannot be assumed to be constant with time. Analysis indicates that it varies in proportion to the square of the flow. Therefore, a Box-Cox transform is used for stabilizing the variance. The appropriate transformation is to take logarithms of the original measurements. The new series has a noise variance which is approximately constant. The method of maximum likelihood is then used to estimate the noise variance for each model. An equivalent log-liklihood function for the flow is written in terms of the observation forecast error; and is evaluated through the use of the Kalman Filter. 5.1 The Estimation Method Because an accurate and reliable estimate of the noise variance is desired, heuristic methods discussed in Chapter 2 are not considered. More sophisticated methods exist, but are traditionally avoided due to extensive computational 54 55 requirements. However, the availability of computer software, has made these methods more attractive. An example is the Bayesian method of estimation. This will not be used in this thesis but is discussed briefly as it has attracted some attention in the recent literature. This method does not estimate the noise variance per se. Actually, it is a discrimination technique. Numerous Kalman models are considered at the start; each one containing a possible value of the noise variance, a2^. It is assumed that one of these Kalman Filters contains the correct noise variance a*2. Posterior probability of a2 ^ being o*2 given the observations, is calculated for each model. After sufficient number of observations, the model with a2- closest to a*2 ' 1 will be assigned the highest posterior probability. Although initial applications by Valdes et al., (1978) show that this method does select the model with a2 ^ closest to the true variance, it is not considered here for the following reasons: 1. Applications of this technique is only in its initial stage, hence its properties still need to be examined. 2. In a study by the previous authors, it 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 original set of models is assumed to encompass the true value of R. 56 This last point presents a limitation in practice. Often, the engineer is not sufficiently familiar with the physical characteristics of a basin to provide the bound values. Hence, the set of a2^'s chosen may not include the true value. Instead, the method of maximum likelihood is used in this chapter to estimate the noise variance. The maximum likelihood method gives estimates which are known to be consistent and asymptotically efficient. No prior knowledge of the system parameters are required here. The difficulty with this method is the evaluation of the likelihood function. In time series applications, errors in successive streamflow measurements are likely to be correlated. Harvey (1981) outlines in his text, how the evaluation of the log likelihood function can be simplified through the use of the Kalman Filter. An example of this estimation procedure is illustrated 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 is 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 will be used in subsequent studies, the noise variance is estimated by the following procedure. 1. Determine the appropriate Box-Cox transformation for the 57 original flow data such that the noise variance of the new series is approximately constant. 2. For the transformed data, obtain the maximum likelihood estimate of the noise variance by using the Kalman Filter. 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 realistic assumption is that the variance of the flow observation error be a function of the flow itself. 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 is assumed. In a real-time forecasting situation, flow prediction is made simpler if the measurement errors are time-invariant. Transforms for stabilizing the variance are available. A Box-Cox transform is used in this application. In general, this class of transforms is appropriate when the variance of the observations can be expressed as a function of its mean. Before the appropriate Box-Cox transformation can be determined, the relation between the variance of the observations and its mean is 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 Filter. The state-space formulation where the ARMAX coefficients are the state variables are considered. The conditional variance of the flow, given the model coefficients, is equal to the variance of the measurement error. This is expressed as the following, where z represents the original flow data measured in m3/s. zfc = H xt 5.1 Var (z, xt) = Var(v ) 5.2 The time subscript is now dropped for notational convenience. It is postulated that the variance of the observation is 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 is 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 local 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 zt+i Mean(z) = 1 s = zfc 5.5 2s+1 Var(z) = i=-s Z (z t + i - V* 2s 1 2s s L 3 i = -s 't + i - z2.(2s+1) 5.6 For s=2, the averaging is with respect to 5 data points; hence this is termed a 5-pt. average. The smaller "s" is, the more local 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 little differences in the values obtained for k. It is decided that a 5-pt. average will be used to obtain local estimates of the mean and variance of the flow. The time period is from April 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 is 1.96. 5.3.3 Conclusions The averaged k value of 1.96 indicates that the variance of the measurement error is approximately proportional to the square of the flow (eqn. 5.3). In other words, the standard deviation of the noise term is directly proportional to the flow. A constant variance for the measurement error would result in k being close to zero. 5.3.4 Box-Cox Transforms A class of transforms used for stabilizing the variance is known as the Box-Cox transforms. An appropriate transformation of the raw flow data can be made. The resulting series will have an error term whose variance is 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 first order Taylor Series expansion of the transformed data is: It is desired that Var (g(z)) be approximately constant. The functional dependence </>, determined in section 5.3.3, is z2. 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 log-likelihood function, log £(y1r y2f««. yT) = log £(y_) is: -| log 2n - |log a2 - \ log | V | - (y_-M)'V"1 (y_-u_) 5.10 Differentiating the above expression with respect to a2 and equating it to 0 yields the maximum likelihood estimate of a2. However, the evaluation of the above expression is time consuming due to |vj; and V"1. This is because V is a non-diagonal matrix as the individual errors are serially correlated. Schweppe (1965) showed that an equivalent expression for log £(y_) can be written in terms of the innovations, v^. It is defined as the observation forecast error, yt-yfc. These innovations are independently normally distributed with mean 0, and variance a2t^, where ffc is calculated by the Kalman algorithm. Hence, £ is a set of independent errors drawn from a multivariate normal, MVN(0, a2D). D is a diagonal matrix defined as: 63 f D = The equivalent log-likelihood expression, log £(y_) is: Because is a linear transformation of y^. which is normally distributed, equations 5.10 and 5.11 are equivalent. The forecast error , and its variance f^. are quantities which can be computed by the Kalman Filter at each time step, as part of the Kalman algorithm. This is achieved by recasting the time series model into the state-space framework. Hence, the Kalman Filter is used as a tool for evaluating the likelihood of T dependent observations from a multivariate normal. Differenting the above expression with respect to a2 yields the formula for a2. -|log 211 - | logo-2 - ^ 2 log f -t=1 T t=1 f 5.11 5.12 The quantities, and ffc are given by the Kalman Filter; where ffc = 1 + HPH'. 64 5.4.2 Results The three ARMAX models given in Chapter 4 will 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 is required. Five estimates are obtained for each model. Flow records from 1976-1980 are used, with the observations taken from April 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. Specification of the noise variance for the Kalman algorithm is based on R_„, in subsequent chapters 6 to 8. 65 5.5 Summary The three ARMAX models described in Chapter 4 will be used for flow forecasting in Chapters 6, 7, 8 respectively. For each model, the variance of the measurement error at each time step is unknown. This chapter is concerned with the estimation of noise variance. In real-time forecasting, it is desired that the variance of the error term be time-invariant; although in practice, this is rarely true. Estimation of the noise variance for each model is divided into two parts and the conclusions are: 1. Regressions of log[Var(flow)] on log[Mean(flow)] indicate that the variance is proportional to the square of the flow. By using the Box-Cox transform, it is determined that the raw flow data should be transformed by ln(flow). The noise variance of the new series is approximately uniform. 2. The method of maximum likelihood is used to estimate the noise variances of the transformed data for each model. The log-likelihood expression is written in terms of the innovations and evaluated by using the Kalman Filter. Average of the maximum likelihood estimates for each model is given. 66 Table 5.3 Averaged Maximum Likelihood Estimate of the Noise Variance Model 1 .0022 Model 2 .0137 Model 3 .0024 6. APPLICATION: AR(1) An autoregressive model of order 1, AR(1) is 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 this study is to compare forecasting performances of the three state-space formulations, based on the following three criteria: 1. RMS error of the flow forecast, 2. maximum relative forecast error, and 3. number of times the forecast error is greater than 25% of the true flow. From the study, the following conclusions are drawn. The best forecasting performance is obtained for the scheme which uses the Kalman Filter to estimate the AR coefficient, a. As the estimate is given resursively, flow prediction is made at each time step with the most recent estimate of a. On the other hand, the Kalman Filter can be used to estimate the flow directly. The formulation which models the flow as the state, and the error term as the observation noise results in the worst forecasting performance. Correction to the flow estimates are not adequate as the Kalman gain matrix is too small. Formulating the AR process as the state equation is equivalent to not using the Kalman Filter at all. The forecasting performance can be improved by casting the model 67 68 in state-space form and applying the filter to estimate the AR coefficient recursively. The model used to describe the streamflow phenomenon at Hope i s: qt = aqt-1 + et 6'1 where q = ln(flow) at Hope a = autoregressive coefficient efc = white noise ~ N(0,a2) 6.1 Description of the Forecasting schemes Three schemes are used to forecast flows at Hope and are described below. The streamflow phenomenon is represented by the AR(1) model of eqn. 6.1. 6.1.1 Properties of Scheme 1 The AR(1) model is recast into state-space framework as follows: State eqn: afc = at_ 6.2 Obs. eqn: qfc = q4._1at + 6.3 The AR(1) process is written as the observation equation, with a is the state variable. Hence, the Kalman Filter is used to give recursive estimates of the AR coefficient. 69 Because time-invariant states are considered, there is no noise term associated with a, and thus Q = 0. The 1-step flow prediction is given by §t+1 = ^t^t+1/t ' "here <*t+l/t ^s based on the latest state estimate given by the filter. Under this scheme, the noise variance is required as input to the Kalman algorithm. 6.1.2 Properties of Scheme 2 The state-space representation of the AR(1) model is: State eqn. qfc = aqt-1 + efc 6.4 Obs. eqn. qfc = qfc 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 this formulation, an estimate of a is required prior to forecasting as it is assumed known by the Kalman algorithm. An approximate least squares estimate, aLg is used. This estimate of a is held fixed throughout the forecasting period. Forecasts are given by the prior estimate of the state vector, Qt+1/t = aLS^t' 6.1.3 Properties of Scheme 3 This final scheme also uses the Kalman Filter to give estimates of the flow directly. However, the AR(1) model is "split up" in the state-space formulation: 70 State eqn: qfc = aqt-1 6.6 Obs. eqn: qmt qfc + efc 6.7 qmt represents the measured flow. Although qfc is the state variable, the noise term is modelled as the measurement noise. Again, an estimate of a is required for the algorithm and aLg is used. Forecasts are given by the filter as Q^+1/t = aLgQt« This differs from Scheme 2 in the expression for qt, as will be shown in section 6.5.4. An estimate for the noise variance is 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) "split up" estimates of a are given recursively by the Kalman Filter estimate of a2 required estimate of a required prior to forecasting equivalent to no Kalman Filtering 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 Filter. The forecasting performance can be affected by the input specification of the Kalman algorithm. For this formulation, these inputs are the initial specifications for the AR coefficient and the noise variance; i.e. a0, P0» R. Schemes 2 and 3 use the Kalman Filter to give estimates of the flow directly, as qfc is the state variable. An estimate of a is required prior to forecasting as it is assumed known by the Kalman algorithm. The main disadvantage of these schemes is that a is 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 is given by the prior estimate of the state qt+1yt = aLSqfc. The latest estimate of the flow is: qt = $t/t_. + Kt[qt-qt] 6.8 The expression for the Kalman gain is: Kt = Pt/t-1H' tHPH'+Rl"1 6.9 where PT/T_1 = a^P^ + Q 6.10 Specification of Q affects Kfc, and thus qt. In this application, H=1 and R=0 leads to a constant gain matrix, Kt=1. Therefore, regardless of the value for Q, Kfc is the same even though pt/t_i changes. Because K=1, equation 6.8 is equivalent to: 73 6.11 Hence the 1-step forecast is: ^t+1/t = aLS qt 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 this situation (Kfc=l), Scheme 2 is equivalent to not using the Kalman Filter. The sensitivity of this scheme's performance is studied with respect to specification of a. Under Scheme 3, the Kalman gain is not equal to 1. In this case, the performance of the filter is expected to depend on the specification of q0, P0, R and a. 6.3 Experimental Procedure The forecasting performance of each scheme is investigated. The data used are streamflow measurements at Hope, from April 1 to September 30, for the years 1981 1982 1983. Their performance is measured by three indicators (PI): 1 . PI 2 74 This relative root mean square error (RRMS) expresses the average magnitude of the forecast error as a fraction. It is measured in terms of the actual units; i.e. y is in m3/s. 2. PI2 = maximum of the absolute relative error, 3. PI3 = number of times the absolute flow forecast error error > 25% of the actual flow. All three indicators are such that the smaller their values, the better the performance. To summarize, PI, represents the average error, PI2 denotes the maximum error, and PI3 indicates the frequency of poor forecasts. Schemes 2 and 3 require estimates for the AR coefficient prior to forecasting. The method of least squares is used to obtain , by treating the time series model as a regression model. An approximate least squares estimate (LSE) for a is the sample correlation coefficient between qt and Qt_i« Taking the dependent variable as qfc, and the independent variable as Qt_ir five regressions are done and the sample correlation coefficient is obtained in each case. Data used are streamflow records at Hope from April 1 to September 30 for the years 1976-1980. The average of the five estimates is used as a for forecasting. 75 Schemes 1 and 3 require the noise variance to be known. Similarly, this has to be estimated prior to forecasting and the method of maximum likelihood is used in this case. The same flow records as above are used. Assuming normality for the observations, the approximate maximum likelihood 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 coefficient The LSE of a obtained for the years 1976-1980, is presented below. Table 6.2 Least Squares Estimates of the autoreqressive coefficient year 1976 1 977 1 978 1 979 1 980 average It is 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 is rounded to .002 for use in the Kalman algorithm. 6.4.3 Performance Indicators for Scheme 1 The filter's performance is tested with respect to a0, P0, and R. The range of parameter values used are given below, followed by the results. 77 Table 6.4 Range of values for input parameters Parameter lower limit upper limit a0 .1 2 P0 .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 all conditions 4% 5% 4% R=2X10"6 4% 5% 5% Table 6.5b Values of PI2 for Scheme 1 Parameter value 1981 1982 1983 all conditions 17% 20% 14% R=2X10-6 18% 20% 12% 78 Table 6.5c Values of PI3 for Scheme 1 Parameter value 1981 1982 1983 all conditions 0 0 0 R=2X10"6 0 0 0 Results show that the filter is very robust to specifications of a0 and P0 for all three performance indicators. It is noted that in this application, the filter is also insensitive to the specification of R. 6.4.4 Performance Indicators for Scheme 2 This scheme is equivalent to forecasting without the Kalman Filter. It was shown in section 6.2 that forecasts do not depend on the value of the noise variance in this application. The forecasting performance is examined with respect to specification of a. The range of values used are: Table 6.6 Range of values for a Parameter a lower limit .95 upper limit 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 PI2 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 PI3 for Scheme 2 Parameter value 1981 1982 1983 a=.95 182 179 182 a=.99 0 1 0 a=1.0 0 0 o=1.01 0 0 0 a=1.05 179 180 182 All three performance indicators show that the forecasting performance is sensitive to specification 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 April 1 to September 30 of each year. Table 6.7c shows that for a=.95 or 1.05, the forecast error is greater that 25% of the actual flow at almost every time step. 6.4.5 Performance Indicators for Scheme 3 The forecasting performance is 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 limit upper limit a .99 1.01 q0 1 15 P0 .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% all values of q0 and P0 46% 48% 44% 82 Table 6.9b Values of PI2 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% all values of q0 and P0 1 30% 77% 1 09% Table 6.9c Values of PI3 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 all values of q0 and P0 116 130 133 All three Performance Indicators show that this is the worst forecasting scheme as even the smallest value for PI, is approximately 50%. Although the performance is robust to 83 specification of q0, P0, and R; it is very sensitive to the specification 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 filter divergence. 6.5 Discussion of Results 6.5.1 Forecasting Performance of Scheme 1 For Scheme 1, the filter's robustness to the AR coefficient, a is not surprising. It was found in Chapter 3 that specification of x0 and P0 has little effect on the forecasting performance. It is also noted that for this application, the filter is also robust to specification of the noise variance. Hence, this formulation of the Kalman Filter is most attractive in real-time forecasting, because a) the forecasting performance is good; and b) the filter is robust to initial specification of model coefficients. Indeed, the main advantage of this 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 correlation coefficients. However, the Kalman Filter does not have this limitation. Consequently, the Kalman Filter can be used to estimate the AR coefficient for non-stationary systems where a is greater than 1. For this application, the estimates of a given by 84 the Kalman Filter are also close to 1. Hence, they support the regression estimates. 6.5.2 Forecasting Performance of Scheme 2 The AR(1) model is written as the state equation of the Kalman model. For K=1, this is equivalent to forecasting without the Kalman Filter. An estimate of a is required prior to forecasting. Results show the filter's sensitivity to specification of the AR coefficient. Values of a outside a 5% range with respect to aLg results in larger PI values. In practice, a cannot be estimated with an accuracy of a few percent. Even slight seasonal variation can cause the true value of a to wander by this amount. 6.5.3 Forecasting Performance of Scheme 3 Forecasting under Scheme 3 results in the worst performance. Large increases in all Performance Indicators are noted when specification of a is changed by as little as 1% from aLg. Therefore, this is not a practical forecasting procedure, as a cannot be estimated with this degree of accuracy. Poor estimates of a once chosen cannot be corrected as it is held fixed during the forecasting period. Flow predictions are given by: $t+l/t = aLS $t 6*13 q. is the Kalman estimate of the state vector for time t: $t = Vt-1 + Kt[(3mt"qi\] 6.14 85 or $t = qt/t-i ^-V + Ktqmt 6-15 There are two limiting cases of interest: 1. Kt=1 results in flow forecasts given by qt+1yt = aLSqmt. This is equivalent to not using the Kalman Filter (Scheme 2). 2. Kfc=0 results in flow forecasts given by Qt+i/t = aLs qt/t-1* Tnis corresponds to not using the measurement information at time t. This scheme leads to propagation of errors as t increases. In fact, this is the reason why this forecasting scheme yields such poor results. 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 q0 raised to an exponent afc, where q0 can be taken to be the initial flow. For large t, 1.01*" results in forecasts, which deviate from the actual flow by a greater amount than .99*". It is shown below that for this application, Kfc is approximately 0. For H=1, the expression for the Kalman gain reduces to: Kt " Pt/t-1 / [Pt/t-1 + R] 6'16 K t is small if pt/t_i is much less than R. Starting with P0, Pyo = a2LsP° + °- ~ p° • 86 Hence, K, ~ p0/[ p0+R], which gives P, ^ (I-K^PQ. 1. If a large P0 is specified, then K,~1 which results in P^O. Subsequent PFC will all be close to 0. 2. If P0 is small compared to R, then K^O. However, P,~P0 which is still small. Hence after a couple of time steps, pt//t_i approaches 0 regardless of P0. 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, is the AR coefficient. Hence the Kalman gain is applied to a. Since small adjustments to a are adequate for improving the observation forecast, only a small Kalman gain is necessary. If the state variable is the flow (as in Scheme 3) then small Kfc values are not sufficient to correct the flow estimates. Therefore, this forecasting scheme is the worst one considered in this chapter because: 1. No adjustment can be made to a. 2. Flow prediction at each time step is not corrected adequately as K is 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 is the best forecasting procedure. This corresponds to formulating the AR(1) model in the observation equation with a as the state variable. Hence, this allows continual updating of the AR coefficient. The scheme is best both in terms of flow forecasts obtained, and robustness to initial specification of a. Scheme 2 formulates the AR(1) process as the state equation in the Kalman model. This is equivalent to not using the filter at all. Although the performance can be as good as Scheme 2, forecasts are sensitive to the specification of the AR coefficient. Scheme 3 is the worst forecasting procedure. Not only are the flow predictions sensitive to a, the performance is always poor. Under this scheme, a small Kalman gain leads to insufficient 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 Filter does not necessarily yield improved flow predictions. The best forecasting scheme is to use the Kalman Filter to estimate the AR coefficient recursively. This allows flow predictions to be made with the latest estimate, afc. 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 is a problem with using the AR(1) model in a Kalman Filter framework to predict observations more than 1 step ahead. In the Kalman model, Hfc which relates the observation to the state, is considered known at time t. For the 1-step forecast, Ht+1 = Qt« But for the 2-step forecast, Hfc+2 = qt+1« The problem is that tomorrow's flow is unknown. In practice, the 1-step forecast is used as an estimate. However, the variance of this forecast error (<3t + 2 ~ Qt + 2^' should not be calculated using the Kalman algorithm. One alternative is to use a different statistical model, which avoids this problem, to predict the flows 2 days in advance. A transfer function model is 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 variable. This input-output model 88 gives the 2-day ahead forecast at daily intervals, as measurements are received everyday. The model used to describe the flow at Hope is a transfer function model with constant coefficients: qHt = [qTt_2 qNt_2] I + et 7. where q = ln(flow) §_ = vector of coefficients e = white noise~N(0,a2=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, Fig. 1.1. The objective is to determine whether flow forecasts given by this model can be improved by using the Kalman Filter to estimate §_ in real-time. Performance of the forecasting schemes are based on the same indicators used for the AR(1) study. The following conclusions are drawn: 1. Flow forecasts are not improved significantly by the Kalman Filter to be of practical interest. 90 2. However, under the Kalman Filter scheme, forecasting performance is robust to specification 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 coefficient as the state yields improved forecasts. Hence, this is the only state-space formulation considered in this study. 7.1.1 Properties of Scheme 1 No Kalman filtering is employed in the first scheme. The transfer function model assumes that the coefficients are time-invariant. Flow 2 days in advance is given by the 1-step forecast: qHt+2 - [QTt qNtl ILS 7.2 where j3LS is the least squares estimate obtained prior to forecasting. 7.1.2 Properties of Scheme 2 Kalman Filtering is used and the transfer function model is cast into state-space format as follows: 91 = "»'t-r = _»'t-i. 7.3 qHt = [QTt_2 qNt_2] lt + efc 7.4 The coefficients j31, |32, are modelled as time-invariant states. Therefore, the Kalman Filter is used to give estimates of £fc recursively at each time step. The flow 2 days in advance is: q"t+2 - [qTt qNt] lt 7.5 where |?t is based on the latest estimate given by the Kalman Filter. 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 Filter would yield better flow forecasts. This is because £ is estimated recursively in real-time. Consequently, flow predictions are made at each time step with the latest estimate of the coefficients. In addition, it is expected that under Scheme 2, the forecasting performance would be robust to the initial specification, j30. This is not likely to be true for Scheme 1. Thus, performance of Scheme 1 is examined with respect to different values of j3rc:« Scheme 2 is examined with respect 92 to input specification of J30 and R. It is not known how the forecasting performance of this proposed model compares with that of the autoregressive model. This is a question of model identification and is not the main objective of this study or thesis. However, the following reasoning can be given. Stations with large drainage areas such as Hope, have relatively slow system response. This is due to the storage characteristics of large basins. Quimpo (1973) showed that this phenomenon can be represented by an autoregressive process. With this consideration, it would not be surprising to find that the forecasting performance of the AR(1) model is 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 is from April 1 to September 30 for 1981 1982 and 1983. The performance of each scheme is measured by the three indicators as defined in Chapter 6. These are: 1. PI, which represents the average error 2. PI2 which denotes the maximum error 3. PI3 which indicates the frequency of poor forecasts Small PI values are indicative of good forecasting performance. 93 Flow prediction under Scheme 1 requires an estimate for J3 prior to forecasting. The method of least squares is used to obtain £Lg. Flow records used for this estimation are based on five years of data, 1976-1980. The same time period of the year is used; April 1 to September 30. An average of the five least squares estimates (LSE) is used as J3 for forecasting. Real-time forecasting with the Kalman Filter requires an estimate for the noise variance, R. Approximate maximum likelihood estimates (MLE) were obtained in Chapter 5 for the years 1976-1980. Again, the average of these five estimates is used as specification for R. 7.4 Results 7.4.1 Estimate of the model coefficients Regression analysis gives values of p2 close to 0 for all years, 1976-1980. Hence only the least squares estimates for |3, are given below. 94 Table 7.1 Least Squares Estimates of the model coefficient 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 significant 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 this input-output model. Thus, p°AV = 1.0613. 7.4.2 Estimate of the noise variance The five 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 is rounded off to .015 for the Kalman algorithm input 7.4.3 Performance Indicators for Scheme 1 Results of the performance measures using j3AV 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 is compared. 96 This corresponds to flow prediction using aLg without the Kalman Filter. Table 7.4 Performance Indicators for the AR(1) model 1981 1982 1983 PI , 5% 7% 6% PI2 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 level a performance which is not as good as that obtained with the autoregressive model. The sensitivity of Scheme 1 with respect to input specification is examined. Streamflow data from 1982 are used as this resulted in the worst performance. Table 7.5 Sensitivity of the Performance Indicators to the model coefficient 0=1.0613 0=1.04 PI, 17% 20% PI 2 48% 46% PI3 24 37 97 The first value of 0 is 0Ay. The second value of 0 is an arbritrary choice which represents a 2% change from the averaged f3AV» There is an increase in PI, (average error) and PI3 (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 PI 1 PI 2 PI , 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 is obtained for the two schemes. In both cases, the number of times the forecast error is in excess of 25% of the actual flow, (PI3) is minimal. The 1982 data reveal a slightly larger difference in the performance indicators between the two schemes. All three performance indicators are larger for Scheme 2 which uses the Kalman Filter. 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 if the Kalman Filter is used. For this scheme, the sensitivity of its performance with respect to the initial specification of j3, and R is 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 limit upper limit initial £ (-.1.9) (1.1) R .0015 15 Under all input specifications, 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 specification of j3. 99 7.5 Discussion of Results Performance Indicators for the two schemes do not clearly indicate which scheme gives better forecasts. However, results of the sensitivity analyses favor the scheme which uses the Kalman Filter. This is because the forecasting performance under Scheme 2 is robust to initial specification of the model coefficents. It is found that without the Kalman Filter, a marked decrease in performance can occur. This was illustrated with a 2% change in the specification of 0 for the 1982 dataset. A regression performed on the 1982 data found that |3LS = 1 .062. This is very close to 0AV 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 is 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 this model is not significantly improved or degraded by applying the Kalman Filter to recursively estimate J3. 2. However, it is recommended that the Kalman Filter be used because the forecasting performance is robust to specification of This is important in situations where the model coefficients cannot be estimated accurately or if they change slowly through time. 8. APPLICATION: ARMAX MODEL The Kalman Filter has been applied to autoregressive and transfer function models to yield estimates of the model coefficients recursively. These models are used for predicting flows at Hope 1 and 2 days in advance respectively. The advantage of using the Kalman Filter is that the forecasting performance is robust to initial specification of the model coefficients. As these are usually unknown in hydrologic applications, this feature is appreciated by practicing engineers. In addition, use of the Kalman Filter gives improved flow forecasts for the AR(1) model. An ARMAX model for describing the flow at Hope is considered in this chapter: qHt - ciqHt_l + CzqTt-2 + c3SNT-2 + et 8'1 This is an AR(1) with upstream inputs. The autoregressive part of this model represents storage characteristics of the river basin, while the exogeneous inputs represent contributions from upstream stations. In this study the Kalman Filter is applied to the ARMAX model to give estimates of the ARMAX coefficients resursively. Flow predictions 1 and 2 days in advance are made. The performance is therefore, measured for both the 1-step and 2-step flow forecasts. The objective is to compare the 1-step forecasting performance of the ARMAX model to that of 101 1 02 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 is 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 is: Ci = "c r c2 = c2 c3 = c3 H r H T N -. , q t - [q t_! q t_2 q t-2] £t + et The Kalman Filter is used to give estimates of the model coefficients c,, c2, and c3 at each time step. Performance is measured for both the 1-step and 2-step flow forecasts. The 1-step forecast is given by: qHt+1 = [qHt qTt_! qNt_i] £t+l/t 8 The 2-step forecast is given by: q t+2 = fq t+1 q t q tJ S.t+2/t 8'5 where £t+1yt and £t+2/t are ^ase<^ on tne latest state 103 estimate given by the filter. As q is unknown, it is 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 is 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. First of all, only one model needs to be used. Secondly, there is only one noise variance to estimate. The station at Hope is characterized by a large drainage area. Thus, it is expected that the autoregressive part of the combined model will dominate. This means that estimates of c, will be relatively larger than those of c2 and c3. With this reasoning, the 1-step forecasting performance should be similar to that of the AR(1) model. The Kalman algorithm requires specification of c 0, P0, and R. From the studies of Chapters 6 and 7, it has been concluded that forecasting performance is robust to initial specification of the model coefficients. Therefore, the sensitivity of the filter's performance is examined with respect to the specification of the noise variance only. One issue of concern here is that the 2-step prediction involves unknown elements in the transition matrix, Hfc+2 of the Kalman model. For the ARMAX model, Hfc+2 = 104 H T N fQ t+1 Q t Q tl« Tne Kalman algorithm assumes that this is known. As tomorrow's flow at Hope is unknown at time t, its 1-step forecast is used as an estimate. Violating this assumption results in an increase in the variance of the forecast error. Therefore, the 2-step forecasting performance of the filter gives an indication of the robustness of the filter to the problem of unknown H in a practical situation. 8.3 Experimental Procedure The forecasting performance is measured by the three indicators described in Chapter 6. For this ARMAX model, the performance is measured with respect to both the 1-step and 2-step forecast errors: 1. 1-step forecast error: Qt+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. PI2 which denotes the maximum error 3. PI3 which indicates the frequency of poor forecast. The nature of these indicators is such that the smaller the values, the better the performance. The forecasting period is from April 1 to September 30 for 1981 1982 1983. 105 For the Kalman algorithm input, an estimate of the noise variance is required, and the approximate maximum likelihood estimates obtained in Chapter 5 are used. 8.4 Results 8.4.1 Estimate of the noise variance Approximate maximum likelihood estimates (MLE) of a2 were obtained in Chapter 5 for the years 1976-1980. The average of these is used as the specification 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 is rounded to .0025 for the Kalman specification. 106 8.4.2 Performance indicators for the ARMAX model Two sets of performance indicators (PI) are given. The first set refers to the 1-step performance of the Kalman Filter. The second set refers to the 2-step performance. 1-step Performance Indicators Table 8.2a Variation of PI, to noise specification 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 Variation of PI2 to noise specification Value of R 1981 1982 1983 .025 26% 17% 21% .0025 34% 17% 24% .00025 35% 17% 24% 107 Table 8.2c Variation of PI3 to noise specification 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 specification 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 PI2 to noise specification 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 PI3 to noise specification 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 relative RMS error for the 2-step forecasts is twice that of the 1-step forecasts. Tables 8.2(a,b,c) show the relative insensitivity of the 1-step PI to specification of R. Table 8.3b indicates a larger decrease in PI2 for the 2-step forecasts when R is 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 Filter is used for updating the model coefficients. 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% PI2 34% 17% 24% 17% 20% 14% PI3 10 0 0 0 0 Values of PI, and PI3 are approximately the same under both models. However, the ARMAX model yields 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 is compared to that of the transfer function model. The formulation used corresponds to letting the model coefficients 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% PI2 104% 32% 34% 40% 54% 33% PI3 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 PI2 are more varied. Forecasting under the ARMAX model yields a higher PI2 value for the 1981 data; while a lower PI2 value is obtained for 1982 data. The 1983 data results in almost the same PI2 values for both models. For PI3, the ARMAX model yields lower values in all three cases. In fact, the number of times the 2-step forecast error is greater than 25% of the actual flow is at most 2. This compares with 32 for the transfer function model. Hence, the ARMAX model seldom gives poor forecasts. 111 8.5 Discussion of Results 8.5.1 Performance of the ARMAX model Tables 8.2c and 8.3c indicate that this model gives relatively few forecasts which result in forecast error greater than 25% of the actual flow. This is true for both the 1-step and 2-step flow predictions. The 1-step forecasting performance is relatively insensitive to specification of the noise variance. However, specififying R larger results in better 2-step performance indicators. This is not surprising as the true noise variance for the 2-step forecast is 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 all three PI are obtained. If one is only interested in predicting the flow one day in advance, then the AR(1) model is adequate. In fact, it is preferred over the ARMAX model because it gives a parsimonious representation. 8.5.3 Comparison of 2-step performance As noted in section 8.4.4, the frequency of poor forecasts (PI3) is less for the ARMAX model than the transfer function model. PI3 is an important indicator if overall consistency and reliability of the forecasts are 112 criteria in practice. In addition, the average forecast error, represented by PI,, for the ARMAX model is less than or equal to that of the input-output model. Hence, in terms of model identification, there is a slight preference for the ARMAX model for predicting streamflow 2 days in advance. The fact that reasonable performance is obtained for the 2-step forecasts with the ARMAX model leads to the following conclusion. For the Kalman model, calculation of the 2-step forecast requires that q is known. As this is unknown in practice, its estimate given by the Kalman Filter as the 1-step forecast is used. Although this violates the assumption of the Kalman model, it is found in this application that the 2-step forecasting performance of the Kalman Filter is acceptable for engineering purposes. In fact, the performance is better than that of the transfer function model. For the hydrologist, it is important that the flow predictions are given with their associated standard error. In the case of the 1-step forecast, this variance is given by the Kalman Filter as HPH' + R. This assumes that H is a deterministic known quantity. When elements of the transition matrix H are unknown, an estimate is used; as illustrated in the 2-step forecast above. Results of the 2-step PI in this study show that reasonable forecasting performance can still be obtained. However, there still remains the problem of what expression is to be used for the variance of the 2-step forecast error. This is addressed in 1 13 the next chapter. 8.6 Summary If flow predictions more than 1-day in advance are required, the ARMAX can be used. This is because it has a comparable performance as the AR(1) model. In addition, the 2-step forecasting performance is better than that of the input-output model. In hydrologic applications, the assumption of known system matrix, H is often not satisfied when predicting future observations more than one time step ahead. This study indicates that using an estimate, H still results in a reasonable forecasting performance. Moreover, it is found that in this particular example, the 1 and 2-step forecasts are quite insensitive to the specification of R. 9. VARIANCE OF THE FORECAST ERROR In flood management, the engineer is often required to predict streamflows several days ahead. The reliability of these predictions is reflected by their associated standard errors. For the Kalman Filter state-space model, *t = **t-l + *t 9*1 Y_t = Htxt + vt 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, is known. In this chapter, a general expression for the variance of the observation forecast error when both H and x are unknown, is developed. ARMAX models have been used for streamflow modelling in this thesis. The state-space formulation considered here is where the ARMAX coefficients are the state variables. The system matrix, H contains past streamflows. The k-step flow forecast, £t+k is given by Ht+k-t+k" For the AR(1) model, flow predictions more than 1 time step ahead requires knowledge of future flows, hence Ht+^ is unknown. For instance, ^t + 2 = fit + 2x-t + 2 = ^-t+]-t + 2 9,3 In practice, the 1-step prediction is 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 this study is 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 subscripts are dropped for notational convenience. 9.1 Background Harrison & Stevens (1976) suggested that the mean square error (MSE) of the k-step forecast can be obtained from the Kalman Filter. For k=1, the filter gives Cov (y_-£) = HPH' + R 9.4 where R and P are covariance matrices for the observation noise, and state estimates respectively. As pointed out by Priestley in the discussion of the above paper (1976), this is not appropriate if H is unknown. Feldstein (1971) developed a formula for Cov (y-y) when both H and x are unknown. Regression models are considered in that context. H corresponds to a design matrix which contains independent variables, and x is a vector of regression coefficients. The formula given by Feldstein assumes that H and x are independent. This assumption presents some limitation in hydrologic applications. As noted for the AR(1) example, H and x are likely to be correlated as they are both based on past values of y. 9.2 Problem definition The statistical model considered is: 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 find 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 its columns one under the other, resulting in a column vector, h. The forecast £ is given by Hx. No assumption is made on the relation between H and x. The covariance matrix for the forecast error is: Evaluation of Cov (Hx) involves calculating the expected value of squared quantities of H and x. This requires the fourth moment of their joint distribution. Expression for the fourth moment can be obtained in terms of lower order moments, if H and x are assumed to be jointly distributed 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 is through the use of the Kronecker product. A preliminary result is given in the next section before the derivation of the general formula. The result contains the expressions for the covariances of products of normal random variables. It is written in terms of the Kronecker product. 9.3 Covariances of Products of normal random variables The Kronecker product is defined as follows: Given 2 matrices A , B ., then A ® B is a (ms,nt) matrix mxn sxt with submatrices consisting of a^jB. Therefore, the Kronecker product of a n-dimensional vector, is a super long vector of length n2. From Magnus and Neudecker (1979), the following result is used: If a is distributed as Nn(y, V) where "n" is the length of a, then Cov (a ® a) is given by (I + Kn) (V ® V + V ® MM' + MM' ® V) 9.7 K is a n2 by n2 matrix defined such that K vec(A) = n J vec(A'). For the purpose of the thesis, the following derivation is restricted to a scalar model for yfc. The AR(1) model is considered where yfc = Yt_1 + e^.. Y^-i ano- at corresPoncl 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 is 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) is given by the 2nd (or 3rd) element of the above vector product. It is 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 is evaluated below, K, = H "2, 1 Z ! 2 identified as X 2 Z2 2 10 0 0 2 0 0 0 0 0 10 (I + K2) = 0 1 1 0 0 10 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 ! 2Z 1 Z2 2 Z ! 2Z Z2 2 2 2 2 2 2 V ® MM' = H ^ 1 1 HxZ, 1 H2Z 1 2 HxZ j 2 x2Z, 1 HxZ 1 2 X 2Z ! 2 H2Z 2 2 HxZ 2 2 X 2 Z 2 2 H L 1 1 H2Z, 2 HxZ 1 1 HxZ ! 2 H2Z2 2 HxZ 1 2 HxZ 2 2 MM' 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 + K2) gives the expression for Var(Hx): H2Z22 + Z„ (x2 + Z22) + £12 <2Hx + ^2) 9.9 Since, Var (y-y) = Var e + Var (Hx) the variance of the forecast error is: (R + H2Z22) + ZV1 (x2 + Z22) + Z12 (2Hx + Z12) 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 is unknown, but its estimate H is uncorrelated with x, then Z12 = 0. This assumption is valid for regression models where the variable H is truly exogeneous. Hence, Var (y-y) = (R + H2Z22) + Z^U2 + Z22) 9.11 This is the same expression as that given by Feldstein's formula for a scalar regression model. 2. If H is known, then Z12 = Z,, = 0. This assumption is used in the Kalman model for calculating the variance of the 1-step innovation. Equation 9.10 then reduces to Var (y-y) = R + H2Z2 2 9. 1 2 Z22 is the mean square error for the state variable x. It is synonymus with "P" in the Kalman Filter notation. 9.5 Illustration of the variance formula for the AR(1) The AR(1) model is used to illustrate the use of the general formula (eqn. 9.10) for the variance of the forecast error. Application of the Kalman Filter to this model was studied in Chapter 6. The state-space formulation considered here corresponds to letting the AR coefficient 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 typically ranges from 1000 m3/s to 7000 m3/s throughout the year. In practice, the AR coefficient is of the order of 1.0. Characteristic values for each term of eqn. 9.10 is examined below: 1 . R Studies from Chapter 5 indicate that the standard deviation of the measurement error is proportional to the flow itself. In practice, these errors are usually less than 10% (Water Survey of Canada, personal communication). Hence, R is taken to be .Oly2. 2 • 2 2 2 This is given by the Kalman Filter as the MSE of the AR coefficient, a. As noted earlier, a is approximately 1. Estimates of a are obtained at each time step from the Kalman Filter. A conservative estimate of its MSE as a percentage of the true value is about 10%. This results in L2 = .01. 3. 2, , This is the MSE for H. For time t+2, this corresponds to yt+1. The MSE(yt+1) is given by the Kalman Filter as HPH' + R. Using the above values for P and R, 2,, = .02y2. 4. 212 Since fl and x are positively correlated, an estimate for Z12 can be obtained by considering the correlation between fl and x. The correlation is of the order of .1 122 to 1, hence I12 is approximately .OOly to .Oly. Each term of eqn. 9.10 is expressed in terms of the flow magnitude: (.Oly2 + .Oly2) + .02y2(1+.01) + .01y (2y + .Oly) 9.13 Disregarding terms with lower orders of magnitude yields (.Oly2 +.0ly2) + .02y2 + .02y2 9.14 Expression 9.14 indicates that in practice, all terms in the general expression are of the same order of magnitude. Therefore, it is not justified to neglect any part of the general formula of eqn. 9.10. Use of the Kalman algorithm to calculate Var (y-y) is equivalent to using the terms in the first bracket only of eqn. 9.14. Thus, the actual variance is significantly larger if the correct expression is used. 9.6 Conclusions A scalar observation model is considered in this chapter, y = Hx + e. A general formula for the variance of the forecast error is developed when H and x are both unknown and their estimates are correlated with each other. The formula obtained is: (R + H2Z22) + In (x2 + 222) + Z12 <2Hx + I12) 123 The expression consists of three parts: The first, accounts for the measurement and state estimation errors (R + £22)* The second part reflects the uncertainty in the estimate, fl (I,,); and the final part accounts for the correlation between fl and x (£12). Two special cases of the general formula are: 1. Feldstein1s formula (R + Z22H2) + I,! (x2 + Z22) This is used if fl and x are uncorrelated which is appropriate if H is truly an exogeneous variable. 2. Kalman Filter formula (R + Z22H2) This is appropriate for the 1-step forecast because H is a known quantity. The importance of this formula is illustrated with an application of the AR(1) model to streamflow prediction at Hope. The variance equation given by the Kalman Filter is inappropriate in this 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 practical illustration indicates that all three parts of the general variance expression are of the same order of magnitude, and hence cannot be neglected. In any forecasting situation, decisions are often made based on the reliability 124 of the future predictions. As the reliability is reflected in the standard error of forecast, the correct expression developed in this 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 Filter to these types of models for flood forecasting is considered in this thesis. The Kalman Filter is a recursive estimation procedure which gives a linear, minimum variance estimator for the state vector at time t. The estimate is 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 this estimation technique depends on satisfying the assumptions of the Kalman state-space model. Use of the Kalman Filter in aerospace applications is 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 this stochastic phenomenon. Choice of the "best" ARMAX model with the "proper" noise statistics is often an impossible goal to achieve in engineering practice. Nevertheless, the Kalman Filter is used even when the assumptions of the Kalman model are not completely satisfied. The following section summarizes the main contributions of the thesis to the field of Kalman Filtering in streamflow 125 126 forecasting. Subsequent sections give the conclusions for each study. The aim of this thesis is to further the understanding of the Kalman Filter as applied to hydrologic systems. The practicality and the performance of this estimation technique is examined in the context of ARMAX flow models. It should be pointed out that the objective is not to identify 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 is for these situations that the application of the Kalman Filter is considered. In this thesis, it is assumed that whichever ARMAX model is chosen, it is an adequate description of the underlying process. 10.1 Summary of Thesis Contributions General practical problems which arise in the application of the Kalman Filter are described in the literature review of Chapter 2. Several problems which frequently occur in hydrologic applications are investigated in this thesis. There are three main contributions of this research, each summarized in the following paragraphs. Advances in the understanding of the Kalman Filter is made in terms of: 1. filter's sensitivity to input specification, and 127 2. practicality of the filter for different state-space formulations of ARMAX models. Correct expression for the mean square error (MSE) of forecast is derived for the AR(1) model. It is shown that under appropriate assumptions, the expression reduces to the Kalman equation for the variance of the innovations. This third contribution is important when autoregressive models are used to forecast flows several time steps ahead. Under the Kalman model, this process violates the assumption of known system matrix H, and use of the Kalman equation underestimates the variance of the forecast errors. The problem of input specification, for quantities often unknown in practice is first investigated. Results indicate that of the four quantities (x0, P0, Q, R), only the combined specification QR has practical effects on the observation forecasts. The sensitivity study gives an insight as to how the noise covariances can be specified in order to achieve reasonable forecasting performances for the filter. In addition, the worst specification combination is also noted. The practicality of the Kalman Filter is illustrated through three special cases of the ARMAX model. The three models are used to describe flow phenomenon of the Fraser River at Hope; typical of basins whose response characteristics are constant, or change slowly through time. The generality of the state-space approach allows flexibility in model formulation. For forecasting purposes, 128 the studies indicate that the most useful formulation is to write the ARMAX model as the measurement equation with the ARMAX coefficients as the state variables. Hence, the Kalman Filter is used to give recursive estimates of the coefficients such that flow forecasts are always made with the latest state estimates. Two main advantages of this formulation over other possible ones are less data requirements, and robustness of flow forecasts to poor initial knowledge of ARMAX coefficients. The Kalman Filter is a 1-step predictor-filter estimation technique. However, forecasts for several time-steps ahead are required in practice and the filter is often used for making these k-step forecasts. In situations where the system matrix H is unknown, the variance of the forecast error should not be calculated from the Kalman algorithm. A correct expression for this variance is developed for the univariate AR(1) model. This expression has important consequences in practice because management decisions are often based on the reliability of flow predictions which is indicated by their mean square errors. It is also shown that in hydrologic applications, all terms in the derived expression are of comparable magnitudes. Conclusions of each investigation are given in the following sections. 1 29 10.2 Sensitivity Analysis The Kalman algorithm requires input specification for the noise covariances and initial conditions of the state vector. As these are usually unknown to the hydrologist, the performance of the filter with respect to misspecification of these inputs (Q, R, x0, P0) is examined. This problem is studied by formulating an ARMAX model in state-space notation with the model coefficients as the state vector. Streamflow data are generated with chosen noise covariances, Q* and R*. Performance of the filter, based on the observation forecast error, is examined with respect to input specifications. It is found that: 1. For the state-space formulation used, initial specification of the model coefficients are not important. Poor choices of x0 and P0 have little influence on the flow predictions. 2. Of the four input factors (Q, R, x0, P0); only the combined specification of the noise covariances have an important effect on the flow forecasts. 3. If Q and R are unknown, it is 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 this state-space formulation, it is important to estimate R correctly. The filter is indifferent to 1 30 misspecification of Q if R is given correctly. 5. However, the reverse is not true. Even if the true Q is given, the filter performs worse if R is 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 likelihood is used to estimate the noise variance. This method is chosen because it gives consistent and asymptotic efficient estimates. Evaluation of the log likelihood function is facilitated by using the Kalman Filter. This is another application of the Kalman Filter, other than that of forecasting. 10.4 The Autoregressive Model The 1-day ahead forecast for the AR(1) model is 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 is to use the Kalman Filter to obtain updated estimates of the AR coefficient. Thus, a is the state variable in the state-space framework. Flow predictions at each time step are made with the latest state estimate. This formulation results in the best forecasting performance. In addition, the flow forecasts are insensitive to initial specification of a. 131 2. Formulating the AR process as the state equation of the state-space model is equivalent to not using the Kalman Filter at all. An estimate of the AR coefficient is required prior to forecasting. Comparable forecasting performance to the best scheme can be obtained for a particular value of a. However, the forecasting performance is sensitive to the specification of the AR coef f icient. 3. The worst performance is obtained for the scheme which "splits up" the AR(1) process in the state-space formulation. Flow is modelled as the state variable while the error term is modelled as the observation noise. Not only is the forecasting performance the worst, it is particularly sensitive to the choice of a determined prior to forecasting. As part of the Kalman Filter algorithm, the 1-step forecast for the observation, and its 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 is known. In practice, a prediction can still be obtained using the 132 1-step forecast from the Kalman Filter as an estimate. The MSE for the 2-step forecast however, cannot be calculated via the Kalman algorithm. The MSE is important as it reflects the hydrologist's confidence in his forecasts. Because costly decisions may depend on the reliability of the flow predictions, the proper expression for the variance of the forecast error should be used. The first approach to this problem is to consider a different ARMAX model which does not require future flows in order to predict the flow 2-days in advance. This is discussed in the next section. The second approach is to derive the correct expression for the mean square error of forecast. This is summarized in section 10.7. 10.5 Transfer Function Model A transfer function model using upstream flow inputs lagged 2 days behind that of Hope is used. Two forecasting schemes are compared, one without and one with the Kalman Filter. The latter scheme formulates the input-output model with the coefficients as the state variables. Conclusions are: 1. For constant coefficients, use of the Kalman Filter does not improve the flow forecasts. 2. However, forecasts obtained without the filter are sensitive to the initial specification of these coefficients. Without the filter, these have to be 133 estimated from past data, thus making the choice of dataset important. 3. Although the filter does not improve the forecasting performance of this model, flow predictions are robust to poor specifications of the model coefficients. Large amounts of data are not necessary for estimating these coefficients prior to forecasting. 10.6 Combined ARMAX Model Using a different model to forecast k steps ahead for various k's involves too many models as k gets large. Therefore a combined model is considered in order to achieve a parsimonious representation of the process. This combined ARMAX model is a combination of the AR(1) and regression models above. The Kalman Filter is applied to this model to give estimates of the ARMAX coefficients 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 Filter is used as an estimate. The study shows that: 1. The 1-step forecasting performance is comparable to that of the AR(1) model. 2. The 2-step forecasting performance is better than that of the transfer function model. Thus, in terms of indentifying the statistical model for 1 34 predicting streamflows up to 2 days in advance, this combined ARMAX model is adequate. The important result of this study is that reasonable forecasting performance is obtained for the 2-step flow predictions; even though an estimate is used for tomorrow's flow in the system matrix H. Violation of the assumption of known H still results in reasonable flow predictions. 10.7 Variance of the forecast error A formula for the variance of the forecast error when both H and x are unknown, is developed for the AR(1) model. It is shown that the new expression calculates a variance significantly larger than that given by the Kalman Filter. This expression for the AR(1) model is: (R + H2Z22) + Z,, (x2 + Z22) + I12 (2Hx + Z,2) All terms in this equation are of the same magnitude in practice. The first part of this expression corresponds to the case of known H, and is equivalent to the Kalman equation for the variance of the innovations. Z22 is the state error covariance matrix. The second part, accounts for the fact that H is unknown but its estimate H is independent of x. Z,, contains the variances and covariances of the elements in H. The final term acknowledges possible correlation between H and x, as reflected in Z12. 135 10.8 State-space Formulation Results of the thesis indicate that application of the Kalman Filter to ARMAX models can give better forecasts when the model coefficients are formulated as the state variables. This is because the filter updates these coefficients recursively in the light of the observation forecast errors. Although it presents some conceptual difficulties, it is necessary to choose between a problem which can be handled in practice and one which is hard to control. The alternate formulation of allowing the flow variable to be the state, requires estimation of the ARMAX coefficients prior to forecasting. This emphasizes the necessity for abundant and good quality data, as forecasting performance is sensitive to the specification of the state transition matrix, 10.9 Future Directions In this thesis, hydrologic systems which can be described by constant coefficient ARMAX models are considered. These are appriopriate when modelling streamflow phenomenon for basins with large drainage areas. Kalman Filtering is applied to these models and the resulting flow predictions are better than or equal to those obtained without the Kalman updating of the model coefficients. A useful extension to this would be to investigate the forecasting performance of the filter for more complex systems whose characteristics change significantly 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. Allin Cornell. Probability, Statistics, and Decision for Civil Engineers. New York: McGraw-Hill, 1970. Bolzern, P., M. Ferraro, and G. Fronza. "Adaptive Real-Time Forecast of River Flow-Rates from Rainfall 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 Filtering." In Applications of Kalman Filter 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 Mifflin Co., 1982. Canada. Environment Canada, Inland Waters Directorate, Water Resources Branch, Water Survey of Canada. Surface Water Data, British 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 Hall, 1980. Chiu, Chao-Lin, and E. 0. Isu. "Application of Kalman Filter to Open Channel Flow Profile Estimation." In Applications of Kalman Filter 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 Filter 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 Filter Modelling of the Speed River Quality." ASCE-J. of the 137 138 Environmental Engineering Division, vol. 105, No. EE5 (Oct. 1979), 961-78. Coulthard, John, ed. UBC Tape. UBC Computing Centre, March 1981. Cox, D. R., and E. J. Snell. Applied Statistics. London: Chapman and Hall, 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 Filter." The Journal of the Institute of Actuaries, vol. 110, part 1, No. 444 (June 1983), 157-81. Dixon, W. J., chief ed. BMDP Statistical Software. Berkeley: Univ. of California Press, 1983. Duong, Nguyen, Rah-Ming Li, and Y. H. Chen. "Adaptive Control of Lock and Dam Gate Openings Using a Kalman Filter for Real-Time Identification of Upstream-Downstream Stage Relationship." In Applications of Kalman Filter 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, vol. 39, No. 1 (Jan. 1971), 55-60. Fox, Daniel J., and Kenneth E. Guire. MIDAS. 3rd ed. Michigan: Statistical 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 Filter to Hydrologyf 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 Statistical 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 Statistics, vol. 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. Tellagraf. version 4. San Diego,CA: ISSCO, 1981. Jazwinski, A. H. "Adaptive Filtering." Automatica, vol. 5 (1969), 475-85. . Stochastic Processes and Filtering Theory. New York: Academic Press, 1970. Johnson, R. A., and Dean W. Wichern. Applied Multivariate Statistical Analysis. Englewood Cliffs, N.J.: Prentice-Hall, 1982. Kahl, Douglas R., and Johannes Ledolter. "A Recursive Kalman Filter Forecasting Approach." Management Science, vol. 29, No. 11 (Nov. 1983), 1325-333. Kalman, R. E. "A New Approach to Linear Filtering and Prediction Problems." Transactions of the ASME - J. of Basic Engineering, (March 1960), 35-45. "A Retrospective after Twenty Years: From the Pure to the Applied." In Applications of Kalman Filter to Hydrology, Hydraulics, and Water Resources. Ed. Chao-Lin Chiu. Pittsburg: Univ. of Pittsburg, 1978. Kendall, M. Time Series. London: Griffin, 1973. Kitanidas, Peter, and Rafael L. Bras. "A Study of Collinearity and Parameter Stability in Rainfall-Runoff Models: Ridge Regression and Kalman Filtering." In Applications of Kalman Filter 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 Filter 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 Bull., 12 (1), (1976), 88-99. Li, Rah-Ming, Nguyen Duong, and Daryl B. Simons. "Application of Kalman Filter for Prediction of State Discharge Relationship in Rivers." In Applications of Kalman Filter to Hydrology, Hydraulics, and Water Resources. Ed. Chao-Lin Chiu. Pittsburg: Univ. of Pittsburg, 1978. Linsley, Ray K. Jr., 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 Filter to Hydrology, Hydraulics, and Water Resources. Ed. Chao-Lin Chiu. Pittsburg: Univ. of Pittsburg, 1978. Lyon-Lamb, Victoria. UBC Tapeasy. UBC Computing Centre, March 1982. McLaughlin, Dennis. "Potential Applications of Kalman Filtering Concepts to Groundwater Basin Management." In Applications of Kalman Filter 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 Statistics, vol. 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 field studies of urban runoff." Canada - Ontario Agreement Research Program, Proj. 73-3-R, Res. Rep. No. 42, 82pp. Mehra, R. K. "On the Identification of Variances and Adaptive Kalman Filtering." IEEE Transactions on Automatic Control, vol. AC-15, No. 2 (April 1970), 175-84. . "Approaches to Adaptive Filtering." IEEE V 141 Transactions on Automatic Control, (Oct. 1972), 693-98. . "Practical Aspects of Designing Kalman Filters." In Applications of Kalman Filter to Hydrology, Hydraulics, and Water Resources. Ed. Chao-Lin Chiu. Pittsburg: Univ. of Pittsburg, 1978. . "Kalman Filters 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. Eric Wood. Oxford: Pergamon Press, 1980. Meinhold, Richard, and Nozer D. Singpurwalla. "Understanding the Kalman Filter." The American Statistician, vol. 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 Statistics with Applications. North Scituate: Duxbury Press, 1973. Moore, John B. WATFIV: Fortran Programming with the WATFIV Compiler. Reston, Virginia: Reston Publishing Co. Inc., 1 975. Moore, Stephen. "Application of Kalman Filter to Water Quality Studies." In Applications of Kalman Filter 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 Filter Estimtation of Parameters and States of Dynamic Water Quality Model." In Applications of Kalman Filter to Hydrology, Hydraulics, and Water Resources. Ed. Chao-Lin Chiu. Pittsburg: Univ. of Pittsburg, 1978. Nakamura, Masatoshi. "Relationship between Steady State Kalman Filter Gain and Noise Variances." Int. J. Systems Science, vol. 13, No. 10, (1982), 1153-163. Natale, L., and E. Todini. "A Stable Estimator for Linear Models." Water Resources Research, vol. 12, No. 4 (Aug. 1976), 667-76. Nicol, Tom, ed. UBC MATRIX. UBC Computing Centre, Mar. 1982. Nightingale, Jon. An Introduction to TEXTFORM. UBC Computing 142 Centre, April 1984. . Textform Layouts. UBC Computing Centre, April 1984. # TEXTFORM. UBC Computing Centre, April 1 984. Patry, Gilles 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 Bulletin, 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 rainfall runoff model." Canadian Journal of Civil 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 Filter in Rainfall-Runoff Studies." In Applications of Kalman Filter to Hydrology, Hydraulics, and Water Resources. Ed. Chao-Lin Chiu. Pittsburg: Univ. of Pittsburg, 1978. Sage, A. P., and G. W. Husa. "Adaptive Filtering with Unknown Prior Statistics." Proc. Joint Automat. Control, 1969, 760-69. Schweppe, "Evaluation of Likelihood Functions for Gaussian Signals." IEEE Transactions on Information Theory, 1965, 61-70. Sittner, Walter T. "WMO Project on Intercomparison of Conceptual Models used in Hydrological Forecasting." Hydrological Sciences Bulletin, XXI, 1 (March 1976), 203-222. Stevens, John, Tony Buckland, and Bruce Joliffe, eds. UBC Fortran. UBC Computing Centre, Nov. 1981. Szollosi-Nagy, A. "An Adaptive Identification and Prediction Algorithm for the Real-Time Forecasting of Hydrological Time Series." Hydrological Sciences Bull., XXI, No. 1 143 (March 1976), 163-76. , Ezio Todini, and Eric Wood. "A State-Space Model for Real-Time Forecasting of Hydrological Time Series." J. of Hydrological Sciences, vol. 4, No. 1 (1977), 61-76. Todini, E. "Mutually Interactive State-Parameter (MISP) Estimation." In Applications of Kalman Filter to Hydrology, Hydraulics, and Water Resources. Ed. Chao-Lin Chiu. Pittsburg: Univ. of Pittsburg, 1978. , and D. Bouillot. "A Rainfall-Runoff Kalman Filter 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 Filter 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 Filter." In Applications of Kalman Filter to Hydrology, Hydraulics, and Water Resources. Ed. Chao-Lin Chiu. Pittsburg: Univ. of Pittsburg, 1978. Walpole, Ronald E., and Raymond H. Myers. Probability and Statistics for Engineers and Scientists. 2nd ed. New York: MacMillan, 1978. Wood, Eric, introd. Real-Time Forecasting / Control of Water Resource Systems. IIASA Proceedings Series. Oxford: Pergamon Press, 1980. Wood, Eric. "An Application of Kalman Filtering to River Flow Forecasting." In Applications of Kalman Filter to Hydrology, Hydraulics, and Water Resources. Ed. Chao-Lin Chiu. Pittsburg: Univ. of Pittsburg, 1978. , and R. K. Mehra. "Model Identification 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, vol. 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. Control, vol. 25, No. 3 (1977), 457-82. Zablosky, J. Paul. UBC Xerox 9700. UBC Computing Centre, 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-06-23
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 |
Aggregated Source Repository | 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:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0062573/manifest