information available at time t, while the intention of smoothing is to include theinformation made available after time t.Harvey (1981) has shown that, if the observations are normally distributed, and thepresent estimator of the state vector is the most accurate, the predictor and the updatedestimator will also be the most accurate. In the absence of the normality assumption asimilar result holds, but only within the class of estimators and predictors which are linearin the observations.Several chapters in this book have documented in detail state space models and theKalman ﬁlter so we will avoid giving a new presentation.7 Let it sufﬁce to say that, basedon the principle of parsimony, our initial Kalman ﬁlter model was an AR(1) process with aconstant mean, implying that extreme temperatures would show some persistence, but theywould eventually return to their mean level for the period under review. Nevertheless,comparing the imputation accuracy of alternative models that we tried, we concludedthat the best model for the Tavg series of the entire dataset is an ARMA(1,2) process,while for the “Autumn” dataset it is an ARMA(2,1) and for the “November” dataset anARMA(1,1).813.3.6 The neural networks regression modelsHere again, two chapters in this book have documented the use of NNR models forprediction purposes.9 In the circumstances, let us just say that the starting point for theNNR models was the linear correlations between the Tmax and Tmin of the index stationconsidered and the explanatory temperature series. While NNR models endeavour toidentify nonlinearities, linear correlation analysis can give an indication of which variablesshould be included. Variable selection was achieved using a backward stepwise neuralregression procedure: starting with lagged historical values of the dependent variableand observations for correlated weather stations, we progressively reduced the number ofinputs, keeping the network architecture constant. If omitting a variable did not deterioratethe level of explained variance over the previous best model, the pool of explanatoryvariables was updated by getting rid of this input (see also Dunis and Jalilov (2002)).The chosen model was then kept for further tests and improvements.7
See Bentz (2003), Billio and Sartore (2003) and Cassola and Luı́s (2003).The Kalman ﬁlter models were implemented with EViews 4.0. The model for the entire dataset is anARMA(1,2) process giving a system of four equations:8
@signal allavg = c(1) + sv1 + c(2)∗ sv2 + c(3)∗ sv3as the observation equation, plus the three state equations:@state sv1 = c(5)∗ sv1(−1) + [var = exp(c(4))]@state sv2 = sv1(−1)@state sv3 = sv2(−1)9
See Dunis and Williams (2003) and Dunis and Huang (2003).
Filling Analysis for Missing Data
391
Finally, each of the datasets was partitioned into three subsets, using approximately 2/3of the data for training the model, 1/6 for testing and the remaining 1/6 for validation.The intentions of this partition are the control of the error and the reduction of the risk ofoverﬁtting. For instance, the ﬁnal models for the “November” Tmax and Tmin series haveone hidden layer with ﬁve hidden nodes.1013.3.7 Principal component analysisPCA is a standard method for extracting the most signiﬁcant uncorrelated sources ofvariation in a multivariate system. The objective of PCA is to reduce dimensionality,so that only the most signiﬁcant sources of information are used. This approach is veryuseful in highly correlated systems, like weather temperatures from neighbouring stations,because there will be a small number of independent sources of variation and most of themcan be described by just a few principal components. PCA has numerous applications inﬁnancial markets modelling: the main ones concern factor modelling within the ArbitragePricing Theory framework (see Campbell et al. (1997)), robust regression analysis in thepresence of multicollinearity (again a problem often affecting factor modelling) and yieldcurve modelling (see Alexander (2001)), although, in this latter case, other techniquesseem preferable.11 It may also be used to overcome the problem of missing data, as weshall see next.13.3.7.1 Principal components computation12Assume that the data for which the PCA is to be carried out consist of M variables indexedj = 1, 2, . . . , M and N observations on each variable, i = 1, 2, . . . , N , generating anN × M matrix X. The input data must be stationary. Furthermore, these stationary datawill need to be normalised before the analysis, otherwise the ﬁrst principal componentwill be dominated by the input variable with the greatest volatility. Thus, we also assumethat each of the M columns of the stationary data matrix X has mean µ = 0 and varianceσ 2 = 1. This can be achieved by subtracting the sample mean and dividing by the samplestandard deviation for each element xij of matrix X. Consequently we have created amatrix X of standardised mean deviations. We will transform this matrix to a new set ofrandom variables, which are pairwise uncorrelated. Let z1 be the new variable with themaximum variance, then the ﬁrst column vector α1 of M elements is deﬁned as:z1(N × 1)
==
X×(N × M) ×
α1(M × 1)
(13.6)
The new variable z1 (N × 1) is a linear combination of the elements in vector α1 . Theproduct zT1 × z1 is the sum of the squares of z1 elements.13Substituting equation (13.6) for z1 we have:zT1 × z1 = (X × α1 )T × (X × α1 ) = α1T × (XT × X) × α110
The NNR models were implemented using the PREVIA software.See Cassola and Luı́s (2003).12Appendix A gives a brief technical reminder on eigenvectors, eigenvalues and PCA.13T indicates matrix transposition.11
(13.7)
392
Applied Quantitative Methods for Trading and Investment
The unbiased covariance matrix of the data which produced matrix X is:1× XT × XN −1where XT × X is (N − 1) times the covariance matrix. Consequently, the maximum variance is found by selecting as vector α1 the one that maximises α1T × (XT × X) × α1 .In addition vector α1 is normalised by the condition that its length equals 1, that isα1T × α1 = 1. In summary, the conditions that determine z1 are:• max zT1 × z1• α1T × α1 = 1
(vector α1 normalised)
This optimisation problem can be solved analytically with the use of the Lagrange multiplier (λ1 ) and vector differentiation. But it can also be solved with Excel Solver.The second vector z2 = X × α2 is determined by the following conditions:• max zT2 × z2• α2T × α2 = 1• α2T × α1 = 0
(vector α2 normalised)(z2 uncorrelated with z1 )
The last constraint derives fromCov(z1 , z2 ) = 0 ⇒ α1T × (XT × X) × α2 = 0 ⇒ (XT × X × α2 )T × α1 = 0⇒ λ2 × α2T × α1 = 0 ⇒ α2T × α1 = 0
(λ2 = 0)
In highly correlated systems, such as weather temperatures, the ﬁrst two eigenvaluesexplain more than 90% of the total variance. There is therefore no need to retain moreprincipal components.13.3.7.2 PCA and missing temperature dataSuppose that the weather station with missing observations is St 1 , and assume thatSt 2 , . . . , St 5 are the correlated weather stations. In order to use PCA to ﬁll in the missingobservations we execute the following steps.• Step 1 : Perform a PCA on St 1 and St 2 , . . . , St 5 using all the available data on St 1to obtain principal components and factor weights (w11 , w12 , w13 , w14 , w15 ). Theselection of the principal components will depend on how highly correlated the systemis. However, they should always be less than the number of variables (<5).• Step 2 : Perform one more PCA on St 2 , . . . , St 5 using all the available data on thesevariables and the same number of principal components (P1 , . . . , P4 ).• Step 3 : Using the factor weights from step 1 and the principal components from step 2we rebuild a simulated data history on St 1 according to the equation:St ∗1 = w11 P1 + · · · + w14 P4• Step 4 : The ﬁnal step is the calibration of our model. The actual data on St 1 that areavailable should be compared with the ones obtained from step 3 to choose the appropriate mean and standard deviation to “reconstruct” the missing weather station data.1414Appendix B documents the Excel ﬁle “WD PCA.xls” of the accompanying CD-Rom, which shows thecomputation of 10 days of missing Tmax temperature data.
Filling Analysis for Missing Data
393
The accuracy of the imputation will depend on the strength of the correlation withinthe system. As weather temperature data are highly correlated, PCA is intuitively anappropriate methodology for ﬁlling missing observations.For our weather observations, all the datasets are stationary, so PCA can be performed.The correlations between temperature data from neighbouring stations are extremely high,so the proportion of total variation that can be explained by only two eigenvalues isabout 96%. For example, the linear model with just two principal components for the“November” Tmax is:Xmax = 0.452 × P1 + 0.1435 × P2The values obtained by the above equation are standardised and therefore need to bemultiplied by their standard deviation and added to their mean to transform them backinto temperatures. We obtain the values for the “November” Tmin data in the same wayand calculate the consequent Tavg temperatures accordingly.
13.4 EMPIRICAL RESULTS13.4.1 The imputation accuracy measuresWe use ﬁve standard error statistics to gauge the imputation accuracy of the alternativeﬁlling methods retained. These measures are documented in Table 13.1 (for more details,see for instance Makridakis et al. (1983)).Table 13.1 Statistical accuracy measuresAccuracy measure
FormulaT1×|ỹt − yt |Tt=1
Mean absolute error (MAE)
MAE =
Mean absolute percentage error (MAPE)
MAPE =
Root mean squared error (RMSE)
T1(ỹt − yt )2RMSE =  ×Tt=1
T100ỹt − yt×Tytt=1
T1 ×(ỹt − yt )2Tt=1
Theil’s inequality coefﬁcient (Theil-U)
U=
Mean Error (ME)
ME =
TT11 ×(ỹt )2 +  ×(yt )2TTt=1t=1T1×(ỹt − yt )Tt=1
394
Applied Quantitative Methods for Trading and Investment
13.4.2 The imputation resultsIn order to conserve space, only the imputed values of the average temperature dataset(Tavg ) are presented and compared. It is worth noting however that, throughout the calculation procedures and for most of the techniques, the imputations of Tavg missing valuesobtained as the average of the Tmax and Tmin imputations for the same missing valuesare more accurate than those obtained by employing the methodologies directly withthe Tavg series. Daily temperatures are reported in degrees Fahrenheit. The comparisonof the imputation accuracy results for the “November” Tavg series of the PhiladelphiaInternational index station (724080) is presented in Table 13.2.These results show that, on all the criteria adopted, PCA substantially outperformed theother methodologies in all cases of missing observations. This means that PCA may beused to avoid mispricing of weather derivatives and to reduce the error in the settlementprocess. The EM, the Kalman ﬁlter and the NNR methods outperform the naı̈ve approach,but are less accurate than the fallback method. In fact, the fallback methodology appears asthe second most accurate method. Finally, the worst methodology in terms of imputationaccuracy appears to be the DA, whereas the naı̈ve approach is validated as a difﬁcult, andTable 13.2Performance Measure
Imputation accuracy measures for the “November” datasetNaı̈ve
Fallback
PCA
Kalman
NNR
EM
DA
One Missing Month Tavg (Average Performance Measures)Mean Absolute Error:8.022.451.136.132.636.228.48Root Mean Squared Error:11.302.961.377.902.977.9011.58Theil’s Inequality Coefﬁcient:0.120.030.010.080.030.080.12Mean Absolute Percentage Error:15.93%4.99%2.35%12.04%5.49% 12.24% 17.50%Mean Error: −0.98−1.68−0.03−2.73−2.57 −2.45 −1.8820 Missing Days Tavg (Average Performance Measures)Mean Absolute Error:8.132.631.006.132.856.257.45Root Mean Squared Error:11.813.221.278.493.188.4911.01Theil’s Inequality Coefﬁcient:0.120.030.010.090.030.090.11Mean Absolute Percentage Error:14.73%5.14%1.99%10.99%5.71% 11.29% 13.74%Mean Error: −4.53−2.03−0.25−4.58−2.75 −4.15 −4.0010 Missing Days Tavg (Average Performance Measures)Mean Absolute Error:4.551.950.653.202.603.455.30Root Mean Squared Error:5.692.550.884.272.904.266.74Theil’s Inequality Coefﬁcient:0.060.030.010.050.030.050.07Mean Absolute Percentage Error:9.85%4.08%1.41%6.68%5.73% 7.28% 11.56%Mean Error:0.751.550.15−0.70−2.400.150.50One Missing Week Tavg (Average Performance Measures)Mean Absolute Error:5.142.570.713.362.713.714.79Root Mean Squared Error:6.363.021.004.593.094.586.07Theil’s Inequality Coefﬁcient:0.070.030.010.050.030.050.07Mean Absolute Percentage Error:10.93%5.31%1.53%6.69%5.90% 7.54% 9.94%Mean Error:2.432.290.29−2.07−2.43 −0.86 −2.07One Missing Day (Tavg )Mean Absolute Error:3.001.000.501.503.503.007.50Root Mean Squared Error:3.001.000.501.503.503.007.50Theil’s Inequality Coefﬁcient:0.030.010.010.020.040.030.08Mean Absolute Percentage Error:6.90%2.27%1.15%3.45%8.05% 6.90% 17.24%Mean Error:3.001.000.50−1.50−3.053.007.50
Filling Analysis for Missing Data
395
Table 13.3 Deviation of computed from actual HDDsMissing Days
ActualHDDs
FallbackHDDs
PCAHDDs
Kalman F.HDDs
EMHDDs
DAHDDs
NaiveHDDs
NNRHDDs
A Missing Day (11/01/93)A Missing Week (11/01/93–11/07/93)10 Missing Days (11/01/93–11/10/93)20 Missing Days (11/01/93–11/20/93)A Missing Month (11/01/93–11/30/93)
21.50123.50186.50293.00487.50
20.00107.50171.00331.00535.50
21.50122.00185.50297.50488.00
23.00138.00193.50378.50563.50
18.50129.50185.00370.00555.00
14.00138.00181.50367.00538.00
18.50106.50179.00377.50511.00
20.00131.50188.50359.50548.00
Deviation from actual HDDs for a Missing Day:Deviation from actual HDDs for a Missing Week:Deviation from actual HDDs for 10 Missing Days:Deviation from actual HDDs for 20 Missing Days:Deviation form actual HDDs for a Missing Month:
1.51615.53848
01.514.50.5
1.514.5785.576
361.57767.5
7.514.557450.5
3177.584.523.5
1.58266.560.5
easy to “compute”, benchmark to beat.15 For the two other datasets, the “Entire” datasetand the “Autumn” temperature observations, the results also show that PCA outperformedall the other methodologies in all the cases of missing data.16Furthermore, using the imputed values of the missing temperature data obtained fromthese techniques, the HDDs for November 1993 have been calculated and compared withactual HDDs. Table 13.3 shows the deviation between computed and actual HDDs. Thebaseline temperature for the computation of the HDDs is 65◦ F.It is again obvious that PCA outperforms all the other methodologies. The deviation ofthe PCA results from actual HDDs for a missing day is 0, while for 10 missing days it is1 degree (against 15.5 degrees for the fallback method); for one missing month, it is 0.5degree, against 48 degrees for the fallback method! On average, the deviation of PCAbased results from actual HDD data is 1.5 degree, while the second best methodology,the fallback method, deviates on average by 23.8 degrees.
13.5 CONCLUDING REMARKSIn this chapter, we set out to examine and analyse the use of alternative methods whenconfronted with missing data, a common problem when not enough historical data orclean historical data exist.Because of the need for “clean” weather temperature data in the fast growing weatherderivatives markets as pricing methodologies rely heavily on them, we examined theimputation accuracy of different interpolation techniques and ﬁlling methods for missinghistorical temperature records. We would argue that the conclusions based on this analysiswould serve as a good guide to the general problem of dealing with missing values.Overall, for the periods and the data series concerned, the results of PCA outperformedall the other methodologies in all cases of missing observations and consequently in thecalculation of HDDs. The only drawback of PCA compared with the second most accuratemethod, the fallback method, is that PCA requires more correlated weather temperature“clean” data. Nevertheless, if the necessary data are available, PCA should be preferredin replacing missing temperature observations. More generally, PCA provides a veryefﬁcient and simple method for ﬁlling missing data in the presence of a correlated systemof variables. As has been shown, it is also easy to implement, as this can be done in Excel.15The poor performance of the EM and DA methods may be linked to the violation of the assumption that thedata are normally distributed (see Section 13.3.4 above), which is not the case here.16Complete results are available from the authors.
396
Applied Quantitative Methods for Trading and Investment
APPENDIX AA.1 Eigenvalue and eigenvectorLet A be a linear transformation symbolised by a matrix A. If there is a vector X ∈ n = 0such that:A×X=λ×Xfor some scalar λ, then λ is called the eigenvalue of A with corresponding (right) eigenvector X. Eigenvalues are also known as characteristic roots, proper values or latentroots (see Marcus and Minc (1988)).Eigenvalues are given by the solutions of the characteristic equation of the given matrix:(A − λ × I) × X = 0where I is the identity matrix.A.2 An introduction to PCAAssume that the data for which the PCA is to be carried out consist of M variables indexedj = 1, 2, . . . , M and N observations on each variable, i = 1, 2, . . . , N , generating anN × M matrix X. As mentioned in Section 13.3.7.1, we also assume that the input dataare stationary and that each of the M columns of the stationary data matrix X has meanµ = 0 and variance σ 2 = 1, which can be achieved by subtracting the sample mean anddividing by the sample standard deviation for each element xij of matrix X. Consequentlywe have created a matrix X of standardised mean deviations.PCA is founded on an eigenvalue and eigenvector analysis of:V = X × X/Nthe M × M symmetric matrix of correlations among the M variables. The principalcomponents are linear combinations of these columns, where the weights are selectedaccording to the following technique:• The ﬁrst principal component describes the majority of the total variation in matrix X,the second component describes the majority of the remaining variation, and so on.• The principal components are uncorrelated with each other.This can be attained by selecting the weights from the cluster of eigenvectors of thecorrelation matrix. If W is the M × M matrix of eigenvectors of V, we have:V×W=W×where  is the M × M diagonal matrix of eigenvalues of V. If matrix W = (wij ) fori, j = 1, 2, . . . , M, then the kth column of W, indicated wk = (w1k , w2k , . . . , wMk ) , isthe M × 1 eigenvector corresponding to the eigenvalue λk . The kth principal componentof the system is given by:Pk = w1k × X1 + w2k × X2 + · · · + wMk × XM
Filling Analysis for Missing Data
397
where Xj indicates the j th column of matrix X, in other words the standardised historicalinput data of the j th variable in the initial correlated system. In matrix notation, the kthprincipal component of the system is given by:Pk = X × wk
APPENDIX BThe ﬁle WD PCA.xls on the accompanying CD-Rom is an illustration of the PCA method.In the “Data” sheet of the workbook, and for the sake of simplicity, we retain 122temperature records in degrees Celsius from 01/09/93 to 31/12/93 for ﬁve weather stations.A 10-day period of missing data is artiﬁcially created between 01/11/93 and 10/11/93.The common dataset for step 1 of the PCA goes until 31/10/93.In the “Step 1” sheet of the workbook, to obtain principal components and factorweights, we perform a PCA on St 1 , St 2 , St 3 , St 4 and St 5 using all the data availableon St 1 up to the missing data period (we have therefore M = 5 and N = 61). We ﬁrststandardise the common dataset up to 31/10/93 (cells H3:L63) by subtracting the samplemean and dividing by the sample standard deviation for each temperature record.To obtain the factor weights, we create two column vectors a1 and a2 with ﬁve cells,one per weather station (initially for a1 , we set the top cell O5 = 1 and cells O6:O9 = 0).We proceed to compute the ﬁrst eigenvector using Excel matrix algebra (cells S4:S64 andcell P14) and the Excel Solver: in the Excel Solver window (Figure B13.1), we maximisecell P14 by changing cells O5:O9 subject to cell P15 = 1.
Figure B13.1
Excel Solver window for the second principal component
398
Applied Quantitative Methods for Trading and Investment
We compute the second eigenvector exactly in the same way: only this time wemaximise cell P18 by changing cells Q5:Q9 subject to cell P19 = 1 and cell P20 = 0.The total variance explained by the ﬁrst two components is computed in cell Q33: at98.4%, there is obviously no need to look for more principal components. We shall retainthe factor weights of St 1 for further use (cells O5:Q5).In the “Step 2” sheet of the workbook, we perform a second PCA, this time only onSt 2 , St 3 , St 4 and St 5 but using the entire databank for these variables and the sametwo principal components (we have therefore M = 4 and N = 122). We shall retain theprincipal components z1 and z2 for further use (cells Q4:Q125 and S4:S125).In the “Steps 3–4” sheet of the workbook, using the factor weights from step 1 (cellsB7:B8) and the useful part of the principal components from step 2 (cells A13:A73 andC13:C73), we rebuild a simulated data history of St 1 according to the equation:St ∗1 = a1 z1 + a2 z2where a1 and a2 are the factor weights for St 1 from step 1 and z1 and z2 the principalcomponents from step 2.In fact, we obtain a standardised series (cells K2:K62) which must be multiplied bythe original series standard deviation and added to its mean to get the “reconstructed”temperature data (cells M2:M62).The choice of the mean and standard deviation is the ﬁnal and crucial calibration stepof the model: here, we take the standard deviation from step 1 (cell E8); for the mean (cellE7), considering both the time dependency of temperature data and the high correlationbetween weather stations, we decided to take the average of the mean temperature forall stations 15 days before and 15 days after the missing data period with the meantemperature for St 2 , St 3 , St 4 and St 5 during the missing data period.Finally, the “Chart” sheet of the workbook plots our results: as can be seen, both actualand PCA-based estimates of the Tavg temperature data for St 1 during the 10-day periodfrom 01/11/93 to 10/11/93 are very close, a view conﬁrmed by the average difference of0.07◦ C over that period (see cell O63 in the “Steps 3–4” sheet of the workbook).
REFERENCESAlexander, C. (2001), Market Models: A Guide to Financial Data Analysis, John Wiley, Chichester.Bentz, Y. (2003), “Quantitative Equity Investment Management with Time-Varying Factor Sensitivities”, Chapter 7 of this book.Billio, M. and D. Sartore (2003), “Stochastic Volatility Models: A Survey with Applications toOption Pricing and Value at Risk”, Chapter 8 of this book.Campbell, J. Y., A. W. Lo and A. C. MacKinlay (1997), The Econometrics of Financial Markets,Princeton University Press, Princeton, NJ.Cassola, N. and J. B. Luı́s (2003), “Modelling the Term Structure of Interest Rates: An Applicationof Gaussian Afﬁne Models to the German Yield Curve”, Chapter 3 of this book.Corbally, M. and P. Dang (2001), “Risk Products”, in E. Banks (ed.), Weather Risk Management:Market, Products and Applications, Palgrave, New York, pp. 105–134.Dempster, A. P., N. M. Laird and D. B. Rubin (1977), “Maximum Likelihood from IncompleteData via the EM Algorithm”, Journal of the Royal Statistical Society, B39, 1, 1–38.Dunis, C. and X. Huang (2003), “Forecasting and Trading Currency Volatility: An Application ofRecurrent Neural Regression and Model Combination”, Chapter 4 of this book.Dunis, C. and J. Jalilov (2002), “Neural Network Regression and Alternative Forecasting Techniques for Predicting Financial Variables”, Neural Network World, 2, 113–139.
Filling Analysis for Missing Data
399
Dunis, C. and M. Williams (2003), “Applications of Advanced Regression Analysis for Tradingand Investment”, Chapter 1 of this book.Greene, W. H. (2000), Econometric Analysis, 4th edition, Prentice Hall, Englewood Cliffs, NJ.Hamilton, J. (1994), Time Series Analysis, Princeton University Press, Princeton, NJ.Harvey, A. C. (1981), Time Series Models, Philip Allan, Oxford.Harvey, A. C. (1989), Forecasting Structural Time Series Models and the Kalman Filter, CambridgeUniversity Press, Cambridge.Henderson, R. (2001), “Pricing Weather Risk”, in E. Banks (ed.), Weather Risk Management: Market, Products and Applications, Palgrave, New York, pp. 167–199.Makridakis, S., S. C. Wheelwright and V. E. McGee (1983), Forecasting: Methods and Applications, John Wiley & Sons, New York.Marcus, M. and H. Minc (1988), Introduction to Linear Algebra, Dover, New York.PricewaterhouseCoopers (2002), “The Weather Risk Management Industry: Survey Findingsfor April 2001 to March 2002