Forecasts of Tropical Pacific Sea Surface Temperatures by Neural Networks and Support Vector Regression by Silvestre Aguilar-Martinez A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF Master of Science in The Faculty of Graduate Studies (Physics) The University Of British Columbia (Vancouver) December, 2008 © Silvestre Aguilar-Martinez 2008 Abstract Nonlinear and linear regression models were developed to forecast the sea surface tem perature anomalies (SSTA) across the tropical Pacific ocean. The methods used were, Bayesian neural networks (BNN), support vector machines for regression (SVR) and lin ear regression (LR). The predictors of the models were a particular combination of the principal components of sea level pressure (SLP) and sea surface temperature (SST) data at lead times ranging from 3 to 15 months. Two data sets corresponding to the time peri ods, 1950-2005 and 1980-2005 were independently studied. The later data set contained the standardized volume of warm water (WWV) integrated across the tropical Pacific basin as an extra predictor. The seasonal cycle was removed from both data sets prior to their introduction into the regression models. Various graphical and statistical tools were used to compare the performance of these models, in particular the correlation scores. The results reveal that although there appeared to be nonlinear elements of the SSTA response to the previous SLP and SST conditions, this response is of second order. The large amplitude variance represented by the first principal component (PC) of the SST field can be adequately expressed as a linear response. Consistently in both data sets, the advantages in the nonlinear modelling manifested mainly through the second PC of the SST variability. The spacial distribution of correlation skills between the predicted and observed SST fields corresponding to the BNN and SVR models show some nonlinear structures in various parts of the domain. Despite some differences in the temporal and spacial distribution of these correlation skills, there was no decisive advantage between using BNN or SVR. Adding the WWV to the predictor set reduced the marginal advantage of the nonlinear methods over the eastern region of the tropical Pacific, in particular at lead times of 9-15 months. These results generally agree with linear recharge oscillator models, and the hypothesis of a linear relationship between thermocline depth and SST in the eastern region. II Table of Contents Abstract ii Table of Contents iii List of Tables v List of Figures vi 1 Introduction 1 1.1 Basic Concepts 1 1.2 Brief History 4 1.3 Atmosphere-Ocean Dynamics 6 1.4 Atmospheric-Oceanic data 11 1.4.1 Sea Surface Temperature 12 1.4.2 Sea Level Pressure 13 1.4.3 Warm Water Volume 13 2 Linear Multivariate Statistical Analysis 15 2.1 Basic Statistical Concepts in Data Analysis 15 2.1.1 Mean, variance and standard deviation 16 2.1.2 Covariance 17 2.2 Principal Component Analysis 17 2.3 Predictors and Predictand 21 2.4 Interpretation of Principal Components 23 3 Statistical Forecast Models 26 3.1 Linear Regression 26 111 Table of Contents 3.2 Neural Networks 28 3.2.1 Training and Optimization 29 3.2.2 Dimensionality 31 3.3 Support Vector Machines 32 3.3.1 Radial Basis Function 36 3.3.2 Implementation of SVR for SST Forecast 36 4 Forecast Results 38 4.1 The 1950-2005 Data Set 38 4.2 The 1980-2005 Data Set 45 4.2.1 Discussion and Conclusion 56 4.3 Further Research 64 Bibliography 65 Appendices A Time Series for El Niño 4,3.4, 3, and 1+2 68 iv List of Tables 4.1 The selected number of hidden nodes is listed for every validation segment. The leftmost entries in the list indicate recent periods, with earlier periods back to 1950 to the right 40 4.2 Hidden nodes results 1980-2005 49 4.3 Values of the hyperparameter set {y, c, e}. The (*) denotes the val ues after WWV was added to the predictor set. The large value of the hyperparameters f(*) reflects the use of nonstandardized pre dictands in the models. The relative large value of € is due to the large amplitude of the PCs (Figure 2.2) 51 v List of Figures 1.1 Niño regions 2 1.2 Recharge Oscillator diagram 9 2.1 The First 5 Modes of SST 24 2.2 The First 5 PCs of SST 25 3.1 feed-forward multilayer perceptron 30 3.2 €-insensitive error norm 33 4.1 1950-2005, correlation scores for the first five PCs of SST . . 39 4.2 Correlation and error for El Niño regions (1950-2005) 43 4.3 BNN: Performance contour map (1950-2005) 44 4.4 SVR: Performance contour map (1950-2005) 46 4.5 1980-2005: correlation scores for the first five PCs of SST . . 47 4.6 Correlation scores for WWV as a single predictor model 48 4.7 1980-2005: Correlation Scores for El Nino 4, 3.4, 3, and 1 + 2. 52 4.8 1980-2005: RMSE Scores for El Nifo 4, 3.4, 3, and 1 + 2 53 4.9 Times series: Niño 3.4, BNN and LR 54 4.10 Times series: Niño 3.4, SVR and LR 55 4.11 BNN: Performance contour map (1950-2005) 57 4.12 SVR: Performance contour map (1950-2005) 58 4.13 BNN: Performance contour map (1950-2005) (+WWV) 59 4.14 SVR: Performance contour map (1950-2005) (+WWV) 60 4.15 LR: Performance contour map 61 A.1 Times series: Niño 4, BNN and LR 69 A.2 Times series: Niflo 4, BNN and LR (+WWV) 70 vi List ofFigures A.3 Times series: A.4 Times series: A.5 Times series: A.6 Times series: A.7 Times series: A.8 Times series: A.9 Times series: A.1O Times series: A. 11 Times series: A.12 Times series: A.13 Times series: A.14 Times series: Niño 4, SVR and LR Niño 4, SVR and LR (+WWV) Nino3.4,BNNandLR Niño 3.4, SVR and LR Niho 3, BNN and LR Niflo 3, BNN and LR (+WWV) Niño3,SVRandLR Niflo 3, SVR and LR (-i-WWV) Niño 1+2, BNN and LR Niño 1 +2, BNN and LR (+WWV) Ninol+2,SVRandLR Nino 1 + 2, SVR and LR (+WWV) 71 72 73 74 75 76 77 78 79 80 81 82 vii Chapter 1 Introduction 1.1 Basic Concepts The El Nino-Southern Oscillation (ENSO) refers to the large scale interaction be tween the tropical Pacific Ocean and the atmosphere at interannual time scales. El Niño/La Nina is the persistent onset of warm/cold surface waters along the tropical central-eastern Pacific Ocean. The Southern Oscillation is the atmospheric com ponent of the ENSO signal, which is characterized by an out-of-phase oscillation of sea level pressure between the western and central-eastern Pacific ocean. Ac cording to the Australian Bureau of Meteorology, the Southern Oscillation Index is defined as the normalized air pressure difference at sea level between the island of Tahiti in the Pacific, and Darwin Australia. SOI= (1.1) Ut1P where LP and op are respectively the mean and standard deviation of the pres sure difference L3.P. The “normal” conditions of the ocean-atmosphere system are defined by the climatology [161. Given the observation record, normal conditions in the tropical Pacific are defined by: • Low pressure over northern Australia • Strong easterly trade winds • Warm waters on the western Pacific and Indian ocean; cold waters on the eastern Pacific off the Peruvian coast • Strong upwelling of cold water along the Peruvian coast Chapter 1. Introduction • Strong storm activity over the warm waters in the west Deviations from this state often results in the development of El Niño, which is defined as a sustained sea surface temperature (SST) anomaly greater than 0.5°C for a period of more than five months across the central tropical Pacific ocean, a region commonly known as Niño 3.4, see Figure 1.1 [1]. The state of ocean- atmosphere system during El Niño is characterized by the following, not exhaustive set of conditions: • Increase in sea level pressure over northern Australia • Weakening of the trade winds • Warm air rises over the Peruvian Coast causing torrential rains • Cessation of upwelling along the Peruvian coast and deepening of the ther mocline 30N 20N IC 1GN NWo1 —‘Zr 2.OS 302OE tOE 1O 150W 120W Figure 1.1: Map of Niño regions La Nina is the complementary phase, when the sea surface temperatures across the central tropical Pacific ocean are unusually low. There is an inherent asymmetry of the ENSO cycle, the cold phase anomalies (La Nina) are typically not as strong as its warm counterpart. This is discussed in more detail in [2]. Depending on the intensity of a El Niño event, its effects can be global and often with catastrophic consequences. For example, the interruption of upwelling 90W 2 Chapter 1. Introduction along the Peruvian coast has a devastating effect on the fishing industry in South America, given that the upwelling of cold, nutrient rich waters is responsible for sustaining the fish populations. Another consequence is an increase in sea level pressure in the western tropical Pacific which can cause the failure of the mon soons in India and the potential for subsequent outbreaks of famine. In fact the famine of 1899 following the failure of the monsoons in India propelled Sir Gilbert Walker, General Director of Observations in India, to investigate the mechanism responsible for the interannual variation of the monsoons [16]. Around the globe El Niño can be associated with a multitude of catastrophes, ranging from wild fires produced by drier than normal conditions, mudslides and floods on noramlly arid regions [16]. These, and other far reaching consequences of El Niflo makes predictions of high practical interest. However, its irregularity in temporal and spacial extensions makes forecasts difficult. A variety of methods have been employed to predict the sea surface temperature anomalies in the tropical Pacific which is the most promi nent feature of the ENSO signal. Physical, statistical and physical-statistical hybrid models have been developed for ENSO forecast [12] [23]. It has been shown that statistical models can attain comparable forecast skill to physical models [20][23]. So far relatively few statistical models have been designed to explore the non linear features of ENSO [23]. An and 3m (2004), suggested that the asymmetry between El Nino/La Nina is a characteristic of the nonlinearities inherent in the ENSO cycle. In this work we have developed a simple linear regression model and two nonlinear statistical models to forecast the SST anomalies over the tropi cal Pacific Ocean. One such nonlinear model is known as Bayesian neural network (BNN) which uses a set of adaptive basis functions to extract the nonlinear relation ships between input data, also known as predictors, and output data (predictand). Due to the general ambiguity associated with neural network training [13], a kernel method known as support vector machine was also used. The kernel methods store the training inputs so that later they can be compared against the set of test inputs. Following this comparison which is characterized by a so called kernel function an output (forecast) can be produced [3]. Chapters 4 and 5 discuss the models and their implementation in more detail. The predictors (inputs) in this case constitute a history of the sea level pressure (SLP) and SST conditions previous to the pre 3 Chapter 1. Introduction diction period. Observations and theoretical models suggest that the warm water volume (WWV) anomaly over the tropical Pacific plays a central role in the ENSO cycle [14] [9] [6]. Therefore, the WWV was included in the models as an extra predictor. 1.2 Brief History In 1889, via a small bulletin article, Dr. Luis Carranza, president of the Lima geographical society, first called attention to a warm counter current flowing north- south along the Peruvian coast. Due to its prominence around Christmas, fishermen from the ports of Paita and Pacasmayo called this counter current “El niño Jesus”, Spanish for the “child Jesus”. Accompanied with particularly intense occurrences of this current was unusually high levels of rainfall over the Peruvian coastal desert and thousands of kilometers inland. These rains brought to the area a rare display of flora and fauna to the otherwise barren desert landscape. The local population began to call these, años de abundancia (years of abundance)[ 161. A few years later, in 1899 the failure of the monsoon rains in India caused a severe famine which prompted Sir Gilbert Thomas Walker to study the interannual variations of the monsoons. His analysis of land and ocean data led to the discovery of an inter- annual seesaw oscillation of atmospheric pressure between the Indian Ocean and the eastern tropical Pacific Ocean [16]. Walker called this fluctuation of air pres sure the “The Southern Oscillation”. Furthermore, He also discovered correlations of the Southern Oscillation with rainfall patterns over India and the tropical Pa cific region [ill. Walker believed that predictability of the Indian monsoon lay on understanding the large scale phenomenon partially revealed by the Southern Oscillation. Inadequate ocean surface data hindered Walker from discovering the effect of the Southern Oscillation on oceanographic variations i.e. sea temperature variations. Walker also documented a correlation between the weakening of trade winds and heavy rainfall in the central equatorial Pacific. It was not until the early 1 940s during a particularly strong aflo de abundancia that the weak trade winds were associated with the warm surface waters that extended from the Peruvian coast far into the central Pacific Ocean. During the 1950s Bjerknes et al. from the University of California, Los Angeles proposed that the exceptionally warm sur 4 Chapter 1. Introduction face waters off the coast of Peru during años de abundancia were a result of large scale changes in the circulation of the entire tropical Pacific Ocean [16]. Bjerknes in 1966 put together the observations obtained after the 1940s El Niflo event to conceive a positive feedback mechanism involving a large scale interac tion between ocean and atmosphere. Simply put, his mechanism uses Walker’s Southern Oscillation to generate a closed circuit of air flow between east and west. Cooled air in the upper atmosphere sinks over the cold surface waters of the east ern tropical Pacific, then zonal pressure gradients forces air to flow westward as trade winds where it collects moisture and rises over the warm western waters as rain clouds. The circuit closes as the flow returns to the east through the upper atmosphere [4]. The mechanism that drives the Walker Circulation depends on heat exchange between the ocean and atmosphere. The temperature gradients over the tropical Pacific determines the air flow, thus interannual variations of the sea surface temperature causes the Southern Oscillation, conversely the Southern Os cillation causes sea surface temperature variations i.e. El Niño/La Nina [4]. This way, Bjerknes was able to summarize many decades of observations into a large- scale ocean-atmosphere interaction, known now as El Nino-Southern Oscillation (ENSO). In subsequent decades, starting in the 1970s, scientist set out to construct var ious dynamical models based on the same fundamental thermodynamical ideas behind Bjerknes’ conceptual model of air-sea interactions. In 1975 Dr. K. Wyrtki from the University of Hawaii published his model on the dynamics of ENSO. He proposed that the built up of warm water on the western tropical Pacific Ocean was caused by “excessively strong southeast trades” [24], two years before the onset of El Niño. This would cause a steep slope of sea level, leading to eastward flow upon relaxation of the wind stress. Since then, many attempts have been made to model the oscillatory nature of ENSO. Many versions of a recharge oscillator model have been published, in which varying mechanisms of the coupled air-sea system pro duce a built up of warm water (charging phase-El Niño). This heat is then dispersed into higher attitudes (discharge phase-La Nina) producing changes in the precipita tion patterns [9] [6]. Other semi-periodic models have tried to incorporate random atmospheric disturbances to explain the spacial and temporal variability of ENSO [161. The most significant source of high frequency atmospheric forcing is the 5 Chapter 1. Introduction Madden-Julian Oscillation (MJO). MJO generates easterward propagating Kelvin waves, then El Niflo conditions in turn enhances MJO [16]. Others, attribute the quasi-periodic and irregular nature of ENSO to various nonlinear processes such as, nonlinear zonal advection, tropical instability waves and even biological physi cal feedback [2]. 1.3 Atmosphere-Ocean Dynamics General circulation models stipulate that the Southern Oscillation (SO) is caused by interannual sea surface temperature variations [16]. Moist air rises over low pressure regions, called convective zones where heavy rainfall is observed [16]. The trade winds blow from high to low pressure regions driving warm surface wa ters to the convective zones which are located where the sea surface temperatures are highest. The tropical Pacific convective zone is located over the western re gion and is responsible for the monsoon rains over India, Indonesia, Malaysia and Northern Australia. This zonal atmospheric circulation, characterized by the air rising over the convective zone and sinking over the regions of high air pressure and low sea surface temperatures is known as the Walker circulation [16]. Changes in location of the warm water has a significant influence on the location of the con vective zones and hence on the large scale rainfall patterns. During El Niño warm surface waters expand eastward, pushing the convective zones to the east. Bjerknes (1969) proposed a mechanism upon which later physical models of El Niño have been built. He stated that a slight relaxation of the trade winds could in turn allow the eastward advection of warm waters. This would further reduce the intensity of the trade winds as the convective zone expands towards the east. The effect is amplified until the trade winds stop, or blow eastward and the warm pool reaches the coast of Peru. Similarly, the effect could start with a shift in the surface waters due to anomalous eastward advection. This positive feedback mechanism rendered the development of various cou pled ocean-atmosphere models [24] [9] [6]. This interaction, according to [16] can be quantitatively described by relating the zonal and meridional atmospheric sur 6 Chapter 1. Introduction face wind stress (TX, r9) to the zonal and meridional (U, V) ocean currents (TX, Tb’) = a(U, V), (1.2) where a is the constant of proportionality. An equation of state is also needed to relate the atmospheric heat source Q, to the temperature of the ocean surface T. Q=/3T, (1.3) with 3 being a constant. As stated in Philander (1990), to relate the atmospheric heat source to ocean dynamics, T must depend on the oceanic variables. The simplest assumed form of this dependence is given by T = sri, (1.4a) Tt+uT—_—rT, (1.4b) where the constant s quantifies the linear dependence between sea surface temper ature, T and j, the depth of the therrnocline. The term uT in (1 .4b) represents mean zonal advection, u is the zonal current speed and T is the temperature gra dient averaged over the meridional direction, and —rT is the damping term with r = constant. Equations (1.4) summarize how changes in the thermocline depth cause a zonal redistribution of warm water and vice-versa [161. Jin (1997) took Bjerknes positive feedback process and the works of Cane, Ze biak and Wyrtki, and convincingly demonstrated from their hypotheses that ENSO is a natural mode of oscillation of the tropical Pacific ocean-atmosphere system [9]. The central idea to most recent models is understanding the role of the ther mocline depth in ENSO dynamics. Assuming mass adjustment of warm water through Sverdrup transport to higher latitudes, the change in thermocline depth in the western equatorial Pacific 1w can be written as: d77WrIlF (1.5) The first term on the right hand side is the damping characterized by the constant r, while the second term, F represents the wind forcing [9]. Furthermore, the 7 Chapter 1. Introduction difference in thermocline depth between east and west is approximately equal to , some value proportional to the zonally integrated wind stress over the equatorial band defined by one Rossby radius of deformation from the equator. Explicitly this can be written as follows, u/E—T1W=T. (1.6) SST is another indicator of the thermocline depth. In the east, water from the subsurface is constantly being upwelled through the climatological trade winds. The weakening of the trade winds and the onset of the SST anomalies produce changes in upwelling and horizontal advection processes. These in turn cause an increase in thermocline depth [9]. These processes can be summarized by the following linear equation. dTE --- —GTE + 7T1E + STE, (1.7) where changes over the averaged central-eastern tropical Pacific SST (TE) are de pendent upon the damping term —CTE and is also proportional to the local ther mocline depth and the averaged zonal wind stress TE over where the SST anomaly takes place [9]. The atmospheric response (, TE) to the SST anomaly TE can be approximated by 1=bTE, (l.8a) TEdTE, (1.8b) where b and dare constants [9]. Substituting (1.8) and (1.6) into (1.5) and (1.7), it can be shown that, dw = —rrlw — GthTE, (1.9a) cITE (l.9b) where a is a constant of proportionality relating and F, and R —c + 7b + öd. Equations (1.9a) and (1.9b) constitute Jin’s linear coupled system between thermocline depth and SST. 8 Chapter 1. Introduction (a) Mature warm phase (c) Mature cold phase (b) Warm to cold transition h (d) Cold to warm transition west east SSTr.,O h Figure 1.2: Recharge oscillator model where (a) and (c) correspond respectively to the warm and cold phases (El Niño/La Nina). The rectangular region represents the tropical Pacific region, with the elliptical region denoting the SST anomaly in the east. The filled arrows describe the wind stress anomaly r, and the unfilled arrows indicate the heat discharge/recharge from the equatorial band. h is the thermocline depth. Panels (b) and (d) represent the transition phases [9]. The physics of this coupled oscillator successfully reproduces some essential features of the ENSO cycle as well as the Bjerknes positive feedback process. Fig ure 1.2 illustrates the essential features of his model. Initially the SST anomaly induces westerly wind forcing which is in phase with the deepening of the thermo dine in the eastern tropical Pacific (Fig 1 .2a). At this point the ocean adjustment west east T +SST west east TO SSTO west east 9 Chapter 1. Introduction processes governed by zonally integrated Sverdrup transport unloads the equatorial band from its heat content resulting in a negative thermocline depth anomaly over the entire tropical Pacific (Fig 1 .2b). The now shallow thermocline permits sea sonal upwelling of cold water in the eastern Pacific. The positive feedback process intensifies the easterly trade winds, leading to anomalous upwelling of cold water in the eastern Pacific, this is known as the cold phase (La Nina) (Fig 1.2c) [9]. This coupled oscillator model is known as the discharge/recharge oscillator model of ENSO. The warm (cold) phase, Figure 1.2a (1.2c) represents El Niño (La Nina) phenomenon. More recently, another oscillator model was proposed by Clarke et al. (2007) where the key variables are the warm water volume (WWV), i.e. the volume of warm water above the 20°C isotherm within the region bounded by [156°E, 80°WJ and [156°S, 5°N], and SST [6]. The physics of Clarke’s model differ slightly from Jin’s recharge oscillator in two aspects. First, according to Clarke’s model the area of maximum ocean-atmosphere coupling takes place in the western-central trop ical Pacific Ocean, contrary to Jin’s where only the eastern SST anomaly drives changes in the thermocline (see equations 1.9). Second, zonally integrated wind stress curl anomalies along the northern and southern boundaries of the box in Clarke’s model produce changes in the WWV by meridional and divergent flow [6]. This mechanism is absent in Jin’s oscillator model. The dynamical variable in Clarke’s model are, the SST anomaly, T and the isotherm depth, i, both aver aged over the western-central region, [156°E, 140°W], [5°S, 5°N]. The coupling of these two variables is given by, (l.lOa) (l.lob) which is Clarke’s couple oscillator model, a and 11 are coupling constants. Mea surements of these constants yields a frequency or w = 27r/(4.25 yr) which is of the right order for ENSO [61. The dynamical models offered by [9] and [6] represent ENSO as a coupled ocean-atmosphere system that features periodic and regular oscillations. The ob 10 Chapter 1. Introduction servations, on the other hand reveal an irregular quasi-periodic cycle with tendency to amplify the positive phase (El Niflo) more than its negative counter part (La Nifla). An and Jin (2004) called attention to the irregularities of ENSO and pro posed a low order nonlinear ENSO model which reproduced the observed asym metries [2]. The heat content in the upper tropical Pacific is given by —a(TE — T0) — — TE), (1.11) ---=-a(Tw-TO)--fj(TW-T) (1.12) where SST is denoted by TE and Tw for the eastern and western region of the tropical Pacific. The other variables include the meridional and vertical velocities (u, w), T0 is the equilibrium temperature, while T’ denotes the temperature of the mixed layer. L and H represent the dimension of the region, L being the width of the basin and H the layer depth [2]. The nonlinearities come in by relations among the wind stress r, temperature (TE, Tw) and velocity (u, w) fields. T —IL(TW—TE) (1.13) u = (1.14) W (1.15) We noticed that this models adds nonlinear coupled dependences between the key variables TE, T, and H. Solutions for this system requires techniques which lie beyond the aims of this project. 1.4 Atmospheric-Oceanic data Following more than half of a century of ocean and atmosphere observations [111, [4] have established that large scale ocean-atmosphere interactions are responsible for the alternating anomalous periods of warm/cold surface waters [16]. The key variables of the ocean-atmosphere system governing the low frequency (interan nual) variations are sea surface temperature (SST) and sea level pressure (SLP). It has been suggested that high frequency atmospheric forcing characterized by 11 Chapter 1. Introduction an index of Madden-Julian Oscillation (MJO) can account for nearly half of the Niño3 .4 variance [141. While the role of high frequency variables such as MJO may be highly significant in accounting for the structural irregularities of ENSO [14], its role was not considered in this work. Instead, two sets of SST and SLP data were used, one covers the period 1950 — 2005, while the other 1980 — 2005. The second period contains warm water volume (WWV) data which is absent from the first set due to its unavailability prior to 1980. 1.4.1 Sea Surface Temperature Sea surface temperature (SST) is perhaps the most relevant feature of El Niño. The SST data used here comes from the NCEP/National Center for Atmospheric Re search (NCAR) reanalysis [10] website, http: / /www. cdc . noaa. gov/ cdc / data . noaa .ersst.html. Monthly average of SST data is available over a global 2° x 2° grid. Version 2 of the extended reconstructed SST (ERSST.v2) used here was available from 1854 to the present [17]. Despite its availability back to 1854, only the periods 1950 — 2005, and 1980 — 2005 were considered. The region of interest is centered in the tropical Pacific (124°E — 70°W, 20°N — 20°S). The data were available in a NetCDF file format; a Matlab program written by graduate student Song Cai of the University of British Columbia, was used to open such data files. Once the data had been opened, smoothing by a three-month running mean was performed, i.e. k+1 Xijk= (1.16) m=k—1 where the variable ijk denotes the smoothed SST data at latitude i, longitude j and time k, and Xjjk is the original SST value. For what followed it was necessary to arrange the gridded variables into columns keeping track of their location on the grid. Each column corresponds to a particu lar month, or more accurately a given three-month average. Two more operations are needed before the data is ready for processing, the first is to eliminate the land data. A given dummy value is assigned to these grid points where SST is not de fined, these points are then removed from the data matrix keeping track of their location, i.e. storing their index. This operation was performed using a matlab file 12 Chapter 1. Introduction specifically written for this purpose. The second operation is the removal of the climatological seasonal cycle associated with each data set. A 57 year seasonal cycle extending through January 1948 to December 2004 was removed from the 1950 — 2005 data set, and a 26 year (January 1980-December 2004) climatology was removed from the 1980 — 2005 data set. The end result is a matrix of SST anomalies where the variables are variations above or below the seasonal mean. Xfl X12 ... Xlm X21 X22 ... X2m (1.17) Xnl Xn2 ... Xnm In this format each xjk represents the SST anomaly at location k on the grid and time 1, the ranges for 1 and k are 1,. . . , n and 1,... , m respectively. Note that the latitude and longitude index (i, j) has been mapped to a single index 1. 1.4.2 Sea Level Pressure Monthly Sea Level Pressure data (SLP) available on a 2.5° x 2.5° global grid was downloaded from NCEP Reanalysis. SLP corresponds to the near surface measurements at (0.995 sigma level) [101. The data used here extends from January 1948 to December 2005 (first set), and January 1980 to December 2005 (second set) over the tropical Pacific ocean (120°E — 70°W, 20°S — 20°N). The data was packaged in the same NetCDF file format as SST, thus was made available for analysis in the Matlab environment. The transformations of running mean, climatology and temporal mean removal described in the previous section were also applied to both sets of SLP data. 1.4.3 Warm Water Volume Warm water volume (WWV) is defined as the volume integral of water above 20°C isotherm between (5°N — 5°S, 120°E — 80°W) [14]. Monthly WWV anoma lies were obtained from the Australian Bureau of Meteorology Research Center (BMRC)OceanAnalysis(http: //www.bom.gov.au/bmrc/ocean/results/ 13 Chapter 1. Introduction climocan . htm) for the period of January 1980 to December 2005. A three month running mean was applied to the data followed by removal of the corre sponding climatological seasonal mean. The index of WWV according to McPhaden et al. (2006) were inferred from the analysis of mooring and expendable bathyther mograph data. 14 Chapter 2 Linear Multivariate Statistical Analysis To undertake a statistical analysis of ENSO, one has to deal with large datasets. A serious study must involve a multitude of variables over a significant time interval. These requirements may in some cases render extraction of meaningful informa tion contained in such datasets a difficult task. A necessary processing of the data must be carried out before it becomes manageable for the purposes of studying ENSO. There is a myriad of available options designed to manage data sets involv ing multiple variables. Among these are: • Regression analysis, linear and nonlinear • Canonical correlation analysis (CCA) • Principal component analysis (PCA) • Nonlinear CCA and PCA Our goal from the outset was to reduce the size of the data structures while retaining important information such as variance. A method for reducing the number of variables while keeping the variance information is principal component analysis (PCA). As described in the next section, this linear method allows for a systematic study of our multivariate SST and SLP data. 2.1 Basic Statistical Concepts in Data Analysis Before one can procede with the development of PCA, it is necessary to lay down some basic concepts and formulae pertaining to statistical analysis. The following 15 Chapter 2. Linear Multivariate Statistical Analysis will not be a rigorous or exhaustive presentation but rather a list of useful and recurrent formulas, for a rigorous treatment, see [211. 2.1.1 Mean, variance and standard deviation The mean of a discrete variable X = {x1, x2, ..., x} as defined by [21] is the expectation value of X given by, E(X) = = xfx(x). (2.1) Here the preferred symbol for the mean will be . The probability function fx (xi) is simply 1/n, assuming all possible values of x are equally probable [21]. Under this simplification, the mean reduces to the simple formula (2.2) Another salient feature of data sets is its dispersion or variance, it measures how ob servations of a variable are distributed about its mean. Mathematically, the variance of a data sample is denoted by the the symbol a2 and it is given by the following formula = (2.3) Sometimes one speaks not of the variable variance but of its square root, known as the standard deviation [21]. Multivariate analysis often requires manipulation of distinct variables into a single regression model. To ensure that differences in units do not produce biased output results, a standardization procedure is necessary. This is accomplished by (2.4) Since a has the same units as x, the above operation also eliminates the units of the input variables [8]. 16 Chapter 2. Linear Multivariate Statistical Analysis 2.1.2 Covariance Clearly the equations presented in the previous section pertain only one dimen sional, or single variable data sets. A multivariate generalization of these concepts is needed for ENSO analysis. In dealing with multiple variables not only is it nec essary to compute the properties of each variable, such as its variance and standard deviation but, we also would like to know how a particular variable changes with respect to changes in another variable. The concept of covariance naturally arises from this need. The covariance between two variables is defined by C(x, y) = (x —X)(Yk — (2.5) Notice that if x = y, the above reduces to the variance formula equation (2.3). These are the essential statistical concepts needed to procede with the development of principal component analysis. 2.2 Principal Component Analysis As explained previously sea surface temperature (SST) is the most important vari able of this study, thus it will be used for the following illustration of principal components analysis. SST over the region of interest can be represented by the set of measurements X(t) = {xi(t),x2(t) ...,x(t)} where t {ti,t, ...,tm}, each value of x represents a measurement of sea surface temperature at a location denoted by i which correspond to a grid point of latitude and longitude. At each location measurements are made a total of m times. This information can be pre sented in the more illustrative form of equation (1.17). At certain times during the development of this project, it was necessary to use an alternative notation to that of equation (1.17), such alternative is, xlj X2j x3= . (2.6) xnj 17 Chapter 2. Linear Multivariate Statistical Analysis Based on this representation, the j index corresponds to a given observation of the set of variables {x3 }, where n is the number of grid points. The covariance matrix for the SST field is = — j)(Xkj — (2.7) rn-i k=1 The advantage of the principal component method is the systematic reduction of the number of variables needed to represent the variance of a data set. This approach follows [8], where, in the language of linear algebra a change of basis is applied to the original data set namely, X in (1.17). The new set of variables denoted by Z is designed to better exhibit the properties of variance. One can think of this geometrically as an axis rotation from {xI,x2,...,x} to {z1,z2,...z}, [8]. The matrix that performs the transformation is e11 e12 ... ei e21 e22 ... e2 E= . . . (2.8) emj em2 ... emn following the notation of [8]. The transformation in matrix and index notation, respectively Z=EX or zik=eijxjk. (2.9) The purpose of this transformation is to find the es such that the variance of z1 is maximized. Recall that the variance of z is m —2 m var(zi)__tZ1) where j =!• (2.10) rn-i m z=1 Substituting (2.9) into yields e1kxk = e1kxk. (2.11) i=1 k=1 k=1 18 Chapter 2. Linear Multivariate Statistical Analysis The variance of zl is var(zi) = (xç — — k=1 1=1 i=1 = e11e1kCkl. (2.12) k=1 1=1 Where we identify Ckl as the covariance matrix defined by (2.7). Using once again the notation of [8] we can write (2.12) in matrix notation var(zi) = eTCei provided e1 = [eu, ..., emiT. (2.13) Recalling that our plan was to chose e1 such that var(zi) is a maximum, however as shown in (2.12) the variance of z is dependent upon the length of ei, to solve this issue the following constraint is imposed [8]. ee1 = 1. (2.14) The technique of choice to solve an optimization problem subject to constraints is the method of Lagrange multipliers. Optimization is performed on the Lagrange function L [8], L = eTCei — — 1). (2.15) One can show that the e1 must be an eigenvector of the covariance matrix and the eigenvalue corresponds to the variance of zi [8]. A = eTCei = var(zi). (2.16) The process is repeated for z2 subject to, ee2 — 1 = 0 with the additional con straint that cov(zi, z2) = 0. It turns out that this condition amounts to the orthog onality relation ee1 = 0. (2.17) The new Lagragian function which incorporates the additional constraints takes the 19 Chapter 2. Linear Multivariate Statistical Analysis following form L = eCe2 — )(ee2 — 1) — 7ee1. (2.18) Finding the stationary points of the Lagrangian function i.e setting the derivatives of L to zero, followed by matrix product with ef on the left leads to an eigenvalue equation for e2 [8]. Therefore, A = eCe2 = var(z2). Previously it was found that the eigenvector e1 which maximizes A = A1, similarly the eigenvector that corresponds to the second largest eigenvalue is A = )‘2. The same procedure is used to find the variances of Z3, Z4, ..., z, which are arranged according to A1 A ... Am 0 and their sum equals the variance of the original variable matrix var(X) = These expressions make evident the primary goal of principal component analysis which is to cast the data matrix into a new set of variables where their characteristic variances are ordered. In some cases only a subset of {z} can account for most of the variance inherent in the set {x }. Our statistical study of ENSO hinges on this special feature of PCA. What remains is to find the coefficients corresponding to the projections in the vector space spanned by the set e X = ETA or, inindexnotation x = a(t)e. (2.19) The projection coefficients a(t) can be easily computed using the orthogonality of condition ee = to yield a(t) = eFx. (2.20) 20 Chapter 2. Linear Multivariate Statistical Analysis To clarify the notation, e11 e21 ... T e12 e22 em2E = : . : = (e1 e2 ... em) and, ei e2fl ... emn ai(ti) ai(t2) . . . a1(t) A = a2(ti) a2(t2) ... a2(tm) (2.21) am(ti) am(t2) ... am(tm) The columns of the matrix ET are the set of eigenvectors of the covariance matrix C. In the literature they are known as principal vectors, spacial patterns, or EOFs (Empirical Orthogonal Functions). The set a (t) are known as the principal com ponents (PCs), or time series. Another important property of PCA is that the PCs { a} are uncorrelated i.e cov(aj,a) = 0, (2.22) [8]. 2.3 Predictors and Predictand As discussed in the previous section, one of the advantages of using PCA is the reduction of the number of variables needed to represent the variance of a multi- variate data set. Explicit use of this property has been made while handling ENSO data. Thus this process defined the variables used in the models. Nevertheless, a presentation of the problem needs to be provided before discussion of the pertinent variables. The solution to the so called regression problem is of primary interest. This problem can be presented via the following equation y = f(x), (2.23) 21 Chapter 2. Linear Multivariate Statistical Analysis where the objective is to find the form of f that better represents the relationship between the variable sets x and y. The first set are called predictors or inputs, the later set, namely y are called predictands (outputs). The relationship f can be linear or nonlinear. The particular choices for f used here are discussed in the next chapter. The predictand, SST in this case, corresponds to the variable whose variability over time one is attempting to forecast. Recall that SST is given by a relatively large data matrix given by equation (1.17), but PCA allows representation of this data by a smaller set of synthetic variables. Therefore, PCA was performed on the SST data matrix (equation 1.17), and the first five PCs served as the predictands in our regression models. Notice that the PCs (a (t)) provide the temporal variations of the variable field, whereas the EOFs (es) account for the spacial variability. Furthermore, since the PCs are uncorrelated, forecasts can be made independently for each one (equation 2.22) [8]. The first five PCA modes of SST account for 59.85%, 10.82%, 7.98%, 3.51%, and 2.85% respectively, of the total variance for the 1950 — 2005 SST data set, and 60.33%, 15.04%, 7.95%, 3.10%, and 2.40% for the 1980 — 2005 data set. The PC time series (at) for 1950-2005 are given in Figure 2.2, while the spacial patterns are given in Figure 2.1. To complete the formulation of the regression problem we need to provide the input variables i.e. x in equation (2.23). The predictors are generated from the PCs of SST and SLP data over a relative short duration prior to the predicted period. In other words, the predictors are the first 9 and 7 PCs of SST and SLP, respectively, with the SLP PCs also supplied at time lags of 3, 6, and 9 months. This procedure incorporates a temporal profile of SST and SLP conditions into the inputs of the regression models [23] [14]. A more specific form of the regression expression appropriate to the forecast models is y(t + ieaci) = f(xi(t),x2(t),x2(t — x2(t —t2),x2(t —t3)), (2.24) where x1 denotes the 9 leading PCs of SST, while X2 denotes the 7 leading PCs of SLP, with lag times t1 = 3 months, t2 = 6 months, and t3 = 9 months. The lead time 1eI as defined in [18] is the time interval between the center of the latest predictor to the center of the predicted period [23]. 22 Chapter 2. Linear Multivariate Statistical Analysis 2.4 Interpretation of Principal Components According to [8], the first three modes of the linear expansion of the SST field by PCA are associated with El NiñofLa Nifla phenomenon. Given that the first three modes account for approximately 75% of the total variance, it is clear that this phenomenon is the dominant feature of the SST variability. Mode 1, (see Figure 2. la) indicates strong variability in the central-eastern equatorial Pacific, increasing from 160°E to the eastern boundary. The peak of the variability occurs at the coast of Peru and extends slightly southward, which is where it was first identified as a warm current flowing south [16]. The time series PCi, Figure 2.2a is the projection coefficient which multiplies every point on the field in EOF1. From Figure 2.2a we can notice that temperatures in the central-eastern tropical Pacific reached their highest levels in 1982/1983 and 1997/1998 El Niño events. Mode 21E0F2 (Figure 2.1b) corresponds to a E-W “standing mode” mentioned in [15] and [22]. This mode reveals the internal structure of SST variability with positive and negative anomalies separated roughly by the 130°W meridian, see Fig. 2. lb. Its time series indicates that PC2 was positive during both El Niflo and La Nina events, and the maximum anomalies also occurred during the two strongest El Niflo events. For mode 3, the peak occurs in the central equatorial Pacific at the 150°W meridian (Figure 2. lc). Shown in Figure 2.2c, the time series of mode 3 illustrates a steady increase in intensity starting at around 1970 [8]. 23 Chapter 2. Linear Multivariate Statistical Analysis ION EQ ‘Os (a) EOF1 59.85% 150€ 180 15(1W 120W 9ôW (b) EOF2 10.82% 150€ 180 150W 120W 90W Figure 2.1: The first five spacial modes (EOFs) for the SST field (1950 to 2005) with positive contours indicated by solid curves, negative contours by dashed curves and zero contours thickened. Contours are multiplied by a factor of 100. 24 -25 c.’J , 20 10 -10 Chapter 2. Linear Multivariate Statistical Analysis Principal Components for SST (modes 1 to 5) 1950 1955 1960 1965 1970 1975 1980 1985 1990 1995 2000 2005 1950 1955 1960 1965 1970 1975 1980 1985 1990 1995 2000 2005 15 10 E5 0 0_ — 0 0 1950 1955 1960 1965 1970 1975 1980 1985 1990 1995 2000 2005 1950 1955 1960 1965 1970 1975 1980 1985 1990 1995 2000 2005 year Figure 2.2: The first five principal components of SST (multiplied 102) for 1950 to 2005. -o 50 0 E 25 0 LC) -o 0 . 5 25 Chapter 3 Statistical Forecast Models The high level of complexity of the ocean-atmosphere system limits a complete analytical development of its evolution. In certain practical cases, we may be more interested in forecasting the SST anomalies than in modeling ocean-atmosphere interactions. According to [20] statistical models yield forecast skills that are com parable to the dynamical models despite their relative simplicity [23]. Many linear statistical models have been used for ENSO forecasting, such as: multivariate lin ear regression, and canonical correlation analysis [23][20]. However, recent evi dence suggests that ENSO may involve nonlinear interactions [2], therefore a non linear statistical treatment is more adequate for such a case. The following sections describe the various statistical models used; the nonlinear methods are Bayesian neural networks (BNN) and support vector machine for regression (SVR). The re sults were compared with multivariate linear regression. 3.1 Linear Regression Linear regression constitutes the simplest nontrivial way to relate two variables, or two sets of variables. Linear regression for a single output variable (SST) can be represented mathematically by y(x,/3) = y(x,) = Xjj+o. (3.1) This expression can be equivalently represented in matrix notation by Y = XT/3 [3]. Recall that by convention, the lower case characters denote a single obser vation of a single variable, while lower case bold characters represent a single 26 Chapter 3. Statistical Forecast Models observation of the full variable set. Bold upper case characters denote a matrix of observations (time series) of the variable set. Generalization to multiple outputs is straightforward and can be done by trading the vector of parameters /3 by a matrix, but since our problem (forecasting SST) involves only a single output we will use the linear model as expressed in equation (3.1). The output y in this regression model is given by a linear combination of the inputs (predictors) {x}. The goal in this case is to choose the set of parameters j3 such that their product with input variables {x } most closely predicts the target value y. To accomplish this, the data set is first divided into training and test sets. The training set constitutes 90% of the time series for every variable; the test set constitutes the remaining 10% of the data. In other words the time series data is split into 10 segments, training is performed using 9 segments and forecasts are made on the remaining segment. The process is repeated until all 10 segments have been forecasted. This procedure is known as cross validation. During every training session the j3s are chosen such that sum of squared residuals is minimized [3]. In order for the following to coincide with the notation used in other texts such as [3], the matrix of inputs X contains n variables along its rows, and each row represents one out of a total of m observations. Then, E = — (3.2) where, = Xjj/3j (3.3) is the model output and y is the target or observed value. Equation 3.2 is mini mized with respect to /3, by 8E/8/3 = 0 [3]. Before preceding it is convenient to transform our variables according to X ‘—* (1, X) and /3 F— (/3, i3) [3]. These transformation allow for equations (3.2) and (3.3) to be cast in a more compact form E = — X/3)2 (3.4) 27 Chapter 3. Statistical Forecast Models The extremum is given by the following = 2(y — Xj3)X1 = 0, i j=O 0. The result can be more easily expressed in matrix notation 2(XTy — XTXI3) 0, /3 = (XTX)_XTy. (3.5) This gives the sets of parameters /3 that minimizes the mean squared error, equation (3.4) [3] for every validation set. 3.2 Neural Networks Attention is now devoted to one of the nonlinear regression models used in the forecast of SST. This model consists of a feed forward two-layer perceptron neu ral network (NN) model with Bayesian regularization [3]. As described in [3], the multilayer perceptron is a particular representation of the general regression equation (3.6) where1(x) corresponds to a set of nonlinear basis functions. In certain large scale applications, the number of basis function can become very large and prohibitively costly in terms of computational resources [3]. An elegant solution is to fix the number of basis function to a relatively small number h, at the cost of allowing the basis functions to depend on some input parameter. This leads to the basic neural 28 Chapter 3. Statistical Forecast Models network structure [31 H = H( wjjx + WjO) where, i = , h. (3.7) Notice that H depends on the parameter set wj and the bias parameter wo. The number or adaptive parameters at the first layer is given by n + 1. Our choice for basis function, following [3] is the hyperbolic tangent function = tanh(wxj + wo). (3.8) These expressions together with the assignments Xj i—* (1, and wj (w0,w) lead to the neural network model used here Yk Wj tanh( wix) + WkO. (3.9) In this case there is only one output variable (SST), therefore only one index is needed to describe the parameters, namely (io,i1). With this model the total number of adaptive parameters to be determined by training are (n + 1)h + 1. The diagram in Figure 3.1 provides a representation of the particular neural network model used in this model. 3.2.1 Training and Optimization Training of the neural network begins with a random choice of the model param eters (w0,w) and (i0,wj) as well as a set number of hidden neurons h. From this choice of parameters an error function such as (3.2) can be computed. How ever, given the flexibility of NN to conform to any continuous function y(x) [3], choosing w, and based on optimization of the error function (3.2) would lead to a problem known as overfitting, which occurs when the network model fits into the noise of the data. To reduce the problem of overfitting, optimization is performed 29 Chapter 3. Statistical Forecast Models xl xO hidden nodes on a regularized error function such as, N E= Yk y1 (3.10) where the second term is a weight penalty or regularization term, and the hyper parameters a and 3 are chosen based on the Bayesian framework, hence the name Bayesian neural network [3]. Equation (3.10) is known as Bayesian regularization, and it was used in the training of the neural network via the MATLAB program “trainbr” from the Neural Network tool box. Finally, the following procedure was used to select the optimum number of hid den neurons (h). As in the case of linear regression the data was partitioned into 10 segment of which 9 were used for training, and one was reserved for testing. For what follows, the labels D9 and Dl are given respectively to the 9 and 1 data seg ments. The D9 set was further partitioned into two sets, one reserved for training and the other for validation. Roughly, 2/3 of the data was used for training, which Hh xri outputs Figure 3.1: General diagram of a feed-forward two layer perceptron. 30 Chapter 3. Statistical Forecast Models by the way was chosen at random in blocks of 12 observations minimum size. This procedure was introduce to further avoid overfitting due to the autocorrelation in the data. The remaining 1/3 of the data was used for validation. For clarity, the labels DT and DV are respectively given to a particular random choice of training and validation data. With the value of h = 1, and a given random pick of DT and DV, the error function (3.10) was optimized a total of 50 times. Since the function (3.10) contains many local minimum, there were potentially 50 different models. The model with lowest mean squared error (MSE) over the training set DT was kept, this error was given the label te. This model was then applied to the corresponding validation set and the MSE over this set was also computed, the name given to this value was Ve. The set (ti, v) characterized this particular choice of DT, DV and h = 1. Next, h was set to 2 and the whole process was repeated without changing DT and DV to obtain (t, v). This process was carried out for /i ranging from 1 to 8, such that in the end there were 8 different sets of pairs (te, ye); the h corresponding to the lowest validation error (ye) is stored as h1. At this point a different DT and DV were randomly chosen from the same D9 and the process was repeated until h2 was determined. This was done a total of 20 times which implies the set s = {h1,h2, ...,h20}. Ideally all the hs should be the same since they come from random picks of the same data set D9. Finally, the mode of set s was chosen as the optimum number of hidden neurons that respresents the data D9 [13]. With this choice of h a total of 50 models were trained using the entire D9 set and tested on Dl. The performance was evaluated in terms of the ensemble average. This entire process was repeated until all 10 segments had been forecasted. 3.2.2 Dimensionality Severe difficulties can arise in problems involving high dimensionality [3]. The problem in a sense results from an exponential increase in the “size” of mathemat ical spaces with an increase in dimensions [3]. In this case, there is a total of 37 time series corresponding to the predictors in the model. The adaptive parameters needed to train such a model must be found in a hypervolume of (37 + 1) >< h + 1 dimensions. The neural network approach used here to forecast SST is therefore 31 Chapter 3. Statistical Forecast Models threatened by the curse of dimensionality. To circumvent this issue, an expedient solution due to [23] was used; the 37 times series were compressed by computing a second PCA, also known as ex tended empirical orthogonal functions (EEOF) analysis, or singular spectrum anal ysis (SSA) [23]. The first 12 SSA PCs were taken to be the final predictors in the model. One Notices that this technique yields a considerable reduction in the dimensions, which in turn yields a reduction of the parameter space volume. The volume is proportional to the number of parameters raised to the power of the num ber of inputs. Therefore, ocy=N25, (3.11) for the special case where N parameters are possible along each dimension, the 37-input parameter space is N25 times larger than the 12-input parameter space. 3.3 Support Vector Machines Support vector machine for regression (SVR) constitutes the other nonliner model used in our forecast of SST. SVR is a particular example of kernel methods [3]. This type of model stores the input vectors during the training process, which are used later during the prediction phase. The method works in the following way. Consider a linear regression model = + b. (3.12) The index i runs over the time series subset of the training data. The set of train ing input vectors x are mapped to a feature space represented by a set of fixed nonlinear basis functions çb (x). Afterwards, linear regression is performed on the feature space characterized by the weight vector w to produce output, [3]. The relative simplicity of expression (3.12) hides the complication involved in evaluat ing a large, possibly infinite set of basis functions k (x). Luckily, there is a way to elude this complication. The SVR method hinges on minimizing a particular regularized error norm defined in terms of a weighted cost function [3]. The cost function contains two 32 where, N flzz1 E((x) — y) = 1 (3.13) (314) Chapter 3. Statistical Forecast Models E6(r) —E +6 r Figure 3.2: c-insensitive error norm, where r = — y. regimes, one called c-insensitive, and linear. For the c-insensitive regime, errors which are smaller than c, where e is a free model parameter, are set to zero, and therefore are not considered. The linear regime assigns a linear cost for errors greater than e, this can be appreciated in Figure 3.2. Mathematically, these condi tion can be summarized by F = C> E(x) — yn) + IJwjJ, Jo, ifI(x)—yI<c — Ill — c, otherwise. The constant C in front of the sum of errors term is an inverse regularization pa rameter, N is the number of observations [3]. The error sum in (3.13) can be expressed only in terms of nonzero error contributions by introducing a new set of variables, called slack variables. These are defined by those points whose predic tion overestimate/underestimate the target value by an amount greater than c. The 33 Chapter 3. Statistical Forecast Models slack variables are denoted by 4 and respectively. An illustration of this idea is given in Figure 3.2. Our regularized error function (3.13) cast in terms of slack variables is F = +i) + IIwII2. (3.15) The variable M is the number of observations for which the input vector x pro duced a target value outside the f-insensitive region. Following Bishop (2006), equation (3.15) must be minimized subject to the constraints 0 and . 0. The problem of optimization subject to constraints lends itself naturally to the use of Lagrange multipliers. Nonzero Lagrange multipliers c, &j, /3, and /3 and a Lagrangian function L are introduced. M M L = F — + — + — yj) (3.16) Substitution of (3.12) into the the above equation followed by setting the deriva tives with respect to w, b, j, and to zero yield the following solution for the model parameters based in terms of the Lagrange multipliers [3]. w = — (3.17) (a — a) = 0, (3.18) (3.19) (3.20) The above conditions are combined to eliminate the variables 3 and 43. The La 34 Chapter 3. Statistical Forecast Models grangian function then takes the form MM L = + &)(j + )T(x.)(x.) i=1 j1 M M — + &j) + + &)yi. (3.21) Clearly, from (3.19), (3.20), and the constraints , & 0 and /3, /3 > 0 we can arrive at the box constraints [3] 0 < a C, (3.22a) 0 < & <C. (3.22b) The set of nonzero multipliers c, and determine the data points x that contribute to prediction. These data points are known as the support vectors and correspond to those points which lie outside the f-insensitive tube [3]. The prediction output is obtained by combining equations (3.17) and (3.12) into (x) = — + b. (3.23) Note that the model output as well as the Lagrangian function (3.21) depend only on the inner product of the basis functions, namely, c5T(xi)(xj) This result, replaces the need to evaluate a large set of complicated nonliner functions with that of defining a proper inner product. In the literature this inner product is known as the kernelfunction. K(x,x3)= (3.24) There is a myriad of kernel functions available, the one used here belongs to a par ticular class traditionally used by beginners called homogeneous kernel functions. These are characterized by a dependence on magnitude of the difference between two input vectors, i.e. K(x,x’) = K(IIx — x’II). 35 Chapter 3. Statistical Forecast Models 3.3.1 Radial Basis Function As discussed in [31, a necessary and sufficient condition for a kernel function K(x, x’), to be valid is that the matrix formed by evaluating K at all possible points in the input space be positive definite. Clearly, this requirement is met if we con sider the linear kernel, K(x, x’) = xTx’. Furthermore, more complicated kernels can be constructed from known valid kernels [3], in particular the exp(K(x, x’)) is valid, provided that K(x, x’) is valid. With this in mind, we consider the following kernel function K(x,x’) = exp(—ilx — x’i12/2u) (3.25) This is known as the Gaussian kernel or radial basis function, for its dependence on lix — x’ . Notice that lix — x’il2 = XTX + 2xTx’ + (x/)Tx/ Expansion of (3.25) in terms of these series of inner products leads to a product of kernels of the type exp(xTx’) which are valid kernels [3]. 3.3.2 Implementation of SVR for SST Forecast The process of implementing SVR was similar to that used in Baysian neural net work and linear regression training. A ten-fold cross validation scheme was de vised to forecast the five leading PCs of SST. The model predictors were the first twelve singular spectrum analysis (SSA) PCs of the combined SLP and SST data set. As before the input data is given by the matrix X, x corresponds to the jth observation of input vectors, where i = 1, 2,... , m. For a given set of observa tions of both input and target values the Lagrangian function (3.21), subject to the box constraints (3.22), was maximized. There are three free hyperparameters in our model which can be adjusted. The value chosen for e determines the number of relevant points used in building the model and thus, determines the number of support vectors set by a1, and d. The inverse of the hyperparameter C controls the amount of weight penalty or regularization, as seen in (3.15). Lastly, the kernel function hyperparameter, denoted by 7 = 1/2u governs the spatial scale of the model, with larger ‘y allowing a finer spatial scale, as ci controls the width of the 36 Chapter 3. Statistical Forecast Models Gaussian kernel in (3.25). The implementation was carried out by a C++ routine interphased with MAT- LAB designed by Chih-Chung Chang and Chih-Jen Lin from the Computer Sci ence department at the National Taiwan University, Taipei. The library of support vector machines (SVM) was available for download at http: / /www. csie. ntu. edu. tw/Thj lin/libsvm. Ten different matlab scripts were written for every cross-validation data segment with the SVM library [5] used to carry out the optimization procedure The training process was conducted by a sevem-fold cross validation procedure for every one of the ten pairs of training and testing sets. A given range of the free parameters e, C, and -y, were probed after which the mean square error (MSE) of each of the seven validation sets evaluated the performance of the model. Finally, the parameters {c, C, -y} that yielded the lowest MSE during validation were chosen for predictions over the test data. This process was con ducted until the entire range had been tested. The search for the hyperparameters was exhaustive, first choosing a grid to identify the region of minimum MSE, then a finer grid for that region until a minimum was identified at the higher resolution level. The results of the SVR method, the Bayesian neural network (BNN) and linear regression (LR) are discussed in the following chapter. 37 Chapter 4 Forecast Results Individual forecasts of the first five PCs of SST were conducted using linear re gression (LR), Bayesian neural networks (BNN) and support vector machines for regression (SVR). Forecasts were made at lead times (Lt1ead) of 3, 6, 9, 12 and 15 months, and the performance of these methods was compared. This analysis was applied to two distinct data records, 1950 — 2005 and 1980 — 2005. The later data set contains an extra predictor, namely, the warm water volume index (WWV). This allowed us to evaluate the importance of WWV on model performance as well as to study the nature of its relationship with the ENSO signal [14][15][6][91. 4.1 The 1950-2005 Data Set For this set, the results of the forecasts of SST support the conclusions reached by Tang et al. (1999). The nonlinear forecast models offer only a marginal advantage over linear regression in the prediction of SST dynamics over the tropical Pacific ocean at seasonal time scales. One possible explanation for the strength for the linear method in comparison to BNN and SVR is that the ENSO signal is basically linear at these time scales. The correlation skills for individual SST PC forecasts are presented in Figure 4.1. The first PC is by far the most significant, accounting for 60% of the to tal variance. Thus, the forecast of PCi is the most important as one would like for the models to successfully predict the time series associated with the largest mode of variability. Recall that PCi has to do with large scale variability over the central-eastern tropical Pacific ocean as it indicates the appearance of warm (or cool) surface waters over Niño 3.4, 3, and 1 + 2 regions (Figure 2.la). As illustrated in Figure 4.1, correlation skills for PCi show no significant difference between the various models at all lead times. The temperature variability in the 38 Chapter 4. Forecast Results r: LR — 0.8 BNNhjii 0.6 r r I lead time (months) Figure 4.1: Cross-validated correlation skill for the 5 leading PCs of SST. The white, black and gray bars denote skill for LR, SVR and BNN, respectively. These scores correspond to the period 1950-2005. eastern-central topical Pacific appears to be dominated by a linear relationship be tween SST and extra-seasonal SLP and SST patterns at 3 to 21 months prior to the predicted month. (a) PCi (b) PC2 (d) PC4 0 C) 3 6 9 12 15 (e) PC5 8 3 6 9 12 15 lead time (months) 39 Chapter 4. Forecast Results The second mode of SST variability associated with PC2 illustrates the east- west contrast in the equatorial Pacific (Figure 2. ib). This second mode was shown in Figure 4. lb to be the most nonlinear aspect of the ENSO signal at short lead times of 3-9 months. Correlation skills between PC2 and the nonlinear forecasts are superior to those of linear regression at these lead times, with BNN show ing the maximum correlation scores. The nonlinearity of PC2 is apparent in the cross-validated selection of the parameter h during model training. Recall that the nonlinearities in the BNN framework are shared between the parameter set {w} and the number of hidden modes h [13]. The only mode for which h = 1 was not selected during the optimization procedure was PC2 at 3, 6, and 9 month lead times. The choice for h in these cases was 2, see Table 4.1. Table 4.1: The selected number of hidden nodes is listed for every validation segment. The leftmost entries in the list indicate recent periods, with earlier periods back to 1950 to the right. PC Leadtime (mon) No. of Hidden Nodes, h Mode 3 1,1,1,1,1,1,1,1,1,1 1 6 1,1,1,1,1,1,1,1,1,1 1 9 1,1,1,1,1,1,1,1,1,1 1 12 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 1 15 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 1 3 2,2,2,2,2,2,2,2,2,2 2 6 2,2,2,2,2,2,2,2,2,2 2 2 9 2,1,2,2,1,2,2,2,2,2 2 12 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 1 15 1,1,1,1,1,1,1,1,1,1 1 3 1,1,1,2,1,1,2,1,1,1 1 6 1,1,1,1,1,1,1,1,1,1 1 3 9 1,1,1,1,1,1,1,1,1,1 1 12 1,1,1,1,1,1,1,1,1,1 1 15 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 Continued on next page 40 Chapter 4. Forecast Results Table 4.1 — continued from previous page PC Leadtime (mon) No. of Hidden Nodes h Mode 3 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 1 6 1,1,1,1,1,1,1,1,1,1 1 4 9 1,1,1,1,2,1,2,1,1,3 1 12 1, 1, 1, 1, 1, 1, 1, 1, 1,3 1 15 1,1,1,1,1,1,1,1,1,1 1 3 2,2,2,2,1,1,1,2,2,2 2 6 2,1,1,1,1,2,1,2,2,1 1 5 9 1,1,1,1,1,1,1,1,1,1 1 12 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 1 15 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 1 As explained in the development of the BNN model, each of the ten segments of the PC time series was forecasted using the remaining 9 segments. The resulting optimum value of h for each of the ten segments is listed in column 3 of table 4.1. The list given in column 3 is also ordered, showing the most recent period on the left the and the more distant past on the right. Only PC2 and PC5 at some lead times selected a model with h higher than the minimum value of h = 1. The PC3 time series corresponds to a mode of variability centered at the Niño 3.4 region [190E-120W, 5N-5S]. This mode accounts for the small out of phase SST oscillation between the Niño 3.4 region and the rest of the domain (Figure 2.1). One feature of this mode is is relatively constant correlation skill from 3 to 12 month lead times. An interesting note about the PC3 time series is its noticeable steady increase over the given time period. The amplitude of this mode was smaller during the 1950s compared to the present with its maximum amplitude occurring during the 1997-98 El Niño episode (Figure 2.2). This observation is in accordance with global warming scenarios and their effect on ENSO, see [7]. These correlation results for PCs 4-5 indicates no advantage of the nonlinear methods over linear regression. The cross validated BNN choice of hidden nodes was 1 for PCs 1, and 3-5 (with the exception of PC5 at 3 month lead time). Recall 41 Chapter 4. Forecast Results that 1 hidden node is the simplest nonlinear model allowed by the neural network method. If h = 1 equation (3.9) reduces to y = tanh(wjx) + o, (4.1) which means there is only one adaptive activation function. Forecasts of the 5 PCs were combined with their corresponding EOFs shown in Figure 2.1 to reconstruct the SST anomalies over the domain. These anomalies were then averaged over the different Niño regions (4, 3.4, 4, and 1 + 2) (Figure 1.1). Correlation scores between the observed and predicted indices was computed for all models as well as the root mean squared error (RMSE). Shown in Figure 4.2, the obtained results are nearly identical for all 3 models. The region with lowest correlation skill corresponds to the Niño 1 + 2 (Fig. 4.2d) which is also the region of largest variability. These results are also very similar to those obtained in [18] where LR, CCA (canonical correlation analysis) and neural networks were compared. Finally, the model performance was evaluated at each grid point over the do main. Observation and predicted time series anomalies were compared via corre lation scores. A contour map of correlation skills for BNN and SVR is provided in Figures 4.3, and 4.4, respectively. These figures also contain the difference in skill between the nonlinear and linear model. As expected, all of our methods show maximum correlation in the central-eastern tropical Pacific at lead times of 3 and 6 months. At longer lead times the correlation skills decrease while show ing maximum values straddling the equator along the 120°W meridian. The right column of panels in Figures 4.3 and 4.4 shows the difference in correlation skills. It was observed, in the case of BNN that at short lead times (3 and 6 months), these differences are positive over most of the domain. In particular, Figure 4.3d resembles spacial mode 2 where the maximum increased correlation regions cor respond roughly to the regions of maximum variability (Figure 2. ib). However, the correlation skill differences are only of the order of 0.01-0.1. Similarly, SVR records a marginally positive difference in the correlation skill over some regions of the domain, particularly over the region centered at 15°S between the 158°E 42 Chapter 4. Forecast Results (a) Nino4 i60E—i5OW, 5N—5S1 (b) Nino4 I ILR 0.8 SVR NN E I FII liii ___________________ (c) Nino3.4 [1 70W—i 20W, 5N—5S] ‘:. III I Hg II i I J I (e) Nino3 [150W—90W, 5N—5S] (f) Nino3 I 3 6 9 nil 5 °. 3 6 9 12 15 (g) Ninoi+2 [90W—80W, ON—i OS] (h) Ninoi÷2 rjiruriI[H1iI lead time (months) Figure 4.2: Correlation scores and root mean square error (RMSE) for the different Niño regions. (d) Nino3.4 3 6 9 12 15 lead time (months) 43 Chapter 4. Forecast Results iON EQ los (a) BNN (Lecid=3 mon) 20S . 150E 180 150W 120W 90W (c) BNN (Lead=6 mon) 1 ON EQ los 150E 180 150W 120W 90W iON EQ los 20S 20N 1 ON EQ los 205 20N 1 ON EQ los (e) BNN (Leod=9 mon) 15OE 180 150W 120W 90W (g) BNN (Lead=12 mon) 150E 180 150W 120W 90W (I) BNN (Lead=15 mon) ThOE 180 150W 12bW 90W (b) BNN — LR 205 150€ 180 150W 120W 90W (d) BNN — LR 9flIJ 150€ 180 150W 120W 9ÔW Figure 4.3: BNN forecast performance and its difference from deviation from the linear models. The contour interval is 0.1 in the right column and 0.02 in the left, with positive (negative) contours indicated by solid (dashed) curves, and the zero contour thickened. --.-/ 1: lOS ION los )\ __ : 44 Chapter 4. Forecast Results and 150°W meridians shown in Figures 4.4d and 4.4f. It also shows a consistent negative correlation difference over the western tropical Pacific for lead times of 3-9 months. Both nonlinear models lose to linear regression at 15-month lead. 4.2 The 1980-2005 Data Set A more in depth analysis was carried out for the shorter data record, namely the 1980-2005 set, due to the availability of a warm water volume index (WWV) from 1980 to the present. The WWV index was included as an extra predictor in all the aforementioned regression models. For comparison purposes two sets of models were developed, one without the WWV index. It was noted that the most intense El Niño events ever recorded occurred during the 1980-2005 period, viz, the 1982- 83 and 1997-98 events (Figure 2.2). The nonlinear terms in equations (1.11) and (1.12) have been found to play a significant role in enhancing the warm cycle dur ing these very strong El Niño events [2]. Hence for the later period of 1980-2005, with more intense ENSO events, the nonlinear methods could be expected to out perform linear ones. As before, the correlation skills for the PCs of SST corresponding to the 1980- 2005 record are given in Figure 4.5. The scores of the linear regression, support vector machines, and Bayesian neural network models are assigned the respective colors, white, black and gray. Two correlation scores are provided for each method. The bar on the left represents the score prior to the inclusion of the WWV index. Before WWV is included as a predictor, BNN provides better correlation scores for the forecast of PCi than SVR, at lead times of 9, 12 and 15 months. SVR shows marginally higher scores than BNN and LR at short lead times (3 and 6 months). Inclusion of the WWV index in the forecast of PCi positively influences the corre lation scores for all models, reducing, or eliminating the differences among them. This change is expected given that the WWV along with SST are the main variables in most theoretical recharge oscillator models [6][9]. Since one is interested in knowing if the SST response to WWV could be gov erned by nonlinear dynamics, the performance of the models was tested using only the WWV index as input. The forecast results for PCi are given in Figure 4.6a, which indicate a linear dependence between SST and WWV, in accordance with 45 Chapter 4. Forecast Results iON EQ ‘Os (a) SVR (Lead=3 man) 205 . — -150E 180 150W 120W 90W (c) SVR (Leod6 mon) 9flN iON EQ los 150E ik 150W l2bW 96W (e) SVR (Leod9 mon) 1 ON EQ los {S 1OE 180 150W 120W 9oW 2ON- (g) SVR (Leadl2 mon) iON EQ ‘Os. 205 20N 1 ON EQ ‘Os 150E 180 i50W 120W 90W (I) SVR (Lead15 mon) 150E 180 150W 120W 96W (b) SVR — LR 20$ 150E 10 150W 120W 96W (d) SVR — LR 9nw .—.... — 20S 1&0E’ ;t;- 120W 96W — (f) SVR — LR 20N .. - iON EQ lOS 150E iái 120W 9ÔW (h) SVR — LR 20N .,.....- 205 150E 180 150W 120W 90W (j) SVR — LR91Th1 . 150E 180 150W 120W 90W Figure 4.4: SVR performance and correlation difference SVR and LR. The contour interval is 0.1 in the right column and 0.D2 in the left. 46 / ION (2.0 7 2 iON iON /-“L2 -. iON -—Ole EQ —02 —608 —0 —000 105 . . . -.-‘ 20S :*.S’ Chapter 4. Forecast Results Figure 4.5: Cross-validated correlation skills for the 5 leading PCs of SST. The white, black and gray bars denote skill for LR, SVR and BNN, respectively. For every method the left(right) bar denotes the score before(after) WWV has been included. These scores correspond to the period 1980-2005. (b) PC2(a) PCi ______ILR 1 I I LR(v) _______ svR ‘ 8 SVR(wwv)BNN BNN(wwv) EJJjj 00 —0. —0.4 (c) PC3 3 6 9 12 15 (d) PC4 0.7 0 C) 8 (e) PC5 3 6 9 12 15 lead time (months) 0 C-) 3 6 9 12 15 lead time (months) 47 Chapter 4. Forecast Results 0 0 (b) PC2 Figure 4.6: Cross-validated correlation skills for PCs 1 and 2 of SST. The white, black and gray bars denote skill for LR, SVR and BNN, respectively. For every method the left bar denotes the score when the model contained the full predictor set (including WWV), whereas the bar on the right denotes the scores when only the WWV was used as a predictor in the models. [15][14]. The correlation scores for PCi are displayed adjacent to the values ob tained using the full set of predictors, with the score bar on the right obtained from the single predictor (WWV) models. For these models the correlation scores pro vided by linear regression were shown to be superior to nonlinear regression at all lead times except for 15 months. Furthermore, the PCi predictions of these models show relatively low scores at very short lead times (3 months). Once again it was noticed that the nonlinearities of ENSO variability come in through PC2. The correlation skills for PC2 at lead times of 3-9 months are supe rior for the nonlinear methods, where SVR attained the best performance (Figure 4.5). The number of hidden nodes selected from the training and validation proce dure for BNN was again 2 for lead times of 3 and 6 months, while all other PCs selected only 1 hidden node (Table 4.2). The hyperparameters chosen during the training of the SVR model are given in Table 4.3. It was discovered that in the case of PC2, inclusion of WWV does not improve the correlation skills (Figure 4.5). The forecast of PC2 using WWV as the single predictor yielded low correlation scores for all models, as indicated in Figure 4.6b. This result suggests that the asymmetry of SST between eastern and western trop (a) PCi I LR I ILR(wv) SVR SVR() BNN BNN(wv) 3 6 9 12 15 3 6 9 12 15 lead time (months) lead time (months) 48 617 ii’i’i’rri’T’i’rici II‘I‘I‘I‘I‘I‘I‘I‘I‘I II’I’I’I’I’I’I’I’I’I6c Ii’I’I’T’i’I’ii’ii9 iT’i’i’ri’T’rii’i II‘I‘I‘I’Ii‘I‘I‘I‘ici II‘I‘1‘1‘1‘1‘T‘I‘I‘111 II’I’f’I’I’I’t’I’t’I6fr II’I’I’I’I’I’I’I’I’I9 II’I’T’t’i’I’I’I’T’I iI’T’z’Ii’I’I’T’z’zci II‘1‘1‘1‘1‘1‘1‘1‘1‘1zI II’I’1’I’I’I’I’I’I’I6£ II’I’I’I’I’I’I’I’I’I9 ii’T’I’rT’I’I’rIi£ II‘I‘I‘I‘I1‘i‘I‘i‘IcI 18‘1‘1‘1‘1‘1‘1‘1‘1‘1 II’T’i’T’I’I’I’I’J’I6Z ZZ’Z.’I’I’Z’I’l’Z’Z’Z9 1l’Z’l’Z’l’Z’I’Z’l’Z£ ‘I‘i‘I1‘111‘11ci II‘I‘I‘I‘I‘I‘I‘I‘I‘I11 II’I’I’I’I’I’I’I’I’I6 II’I’I’I’I’t’I’T’I’I9 II’I’I’I’1’I’I’I’I’I£ POJAJ‘dSPONuppjjooN(uoul)wJpn3j cOO-Og6iSpOUUPP!Hqj S3dsqiijji1CI1UY!U!S:ouP!P)(PU!AMMqi JOUO!snpu!1qM‘c171flIU!UA!as3Jojsios(AM)qdp uuposi3otjjoiupudpuiii!JO1qpozuopiwqoO!JIDId sjjnsajsaIoJj Chapter 4. Forecast Results Following the procedure of the previous section, the predicted five PCs were combined to generate the Niflo 4, 3.5, 3, and 1 + 2 indices. The correlation skills and root mean squared error (RMSE), before and after WWV was introduced into the input set are given in Figures 4.7 and 4.8, for the various Niño regions. Before the WWV index was added BNN provided a better correlation skill than SVR at lead times of 9-15 months for the Niño 3.4 and 4 regions. Incidentally, correlation skills for the Niflo 4 region show the largest observed discrepancies among the dif ferent models, with BNN dominating at longer lead times (9-15 months), and SVR at short lead times (3 and 6 months). Incorporating the WWV tended to reduce the differences in correlation among the 3 methods due to the linear coupling between WWV and SST [6]. For region Niño 3.4 the differences in correlation scores be fore and after the addition of WWV are less pronounced. The superiority of the nonlinear methods is more apparent in the region of largest variability, Niño 1 + 2, over short lead times of 3-9 months. However, the correlation scores for LR show improvements at long lead times for all regions except Niño 4. This is consistent with the forecast results of PCi (Figure 4.6a), where LR provides better correlation scores at 6-15 months lead time. An alternative presentation of the results for the various Nifio regions is given by plotting simultaneously the time series of observed values and the predicted values of LR and one of the nonlinear models. Figures 4.9 and 4.10 illustrate the time series of the observed El Nub 3.4 compared to the predicted time series of the LR, BNN and SVR models (without WWV as a predictor). Notice that at 3 month lead time all of the models successfully predict the major El Niño episodes of 1982/1983, 1986/1987, 1991/1992, and 1997/1998. Beyond 9 months lead, identification of anomalies becomes difficult, however, the BNN prediction and the observed anomaly closely match for El Niño 1982/1983, 1986/1987, and La Nina 1984/1985, 1988/1989. The time series for the other 3 Niño regions and for Niflo 3.4 with WWV added are provided in Appendix A. The SST anomalies were reconstructed by combining the PCs forecasts with their corresponding EOFs. Contour maps of the correlation scores between the predicted and observed anomalies for the nonlinear models are provided in Figures 50 Chapter 4. Forecast Results Hyperparameters PC Leadtime ‘y C (x103) C (x103) 3 0.015 25.00 512 0.011 12.50 60 6 0.041 9.10 416 0.011 25.00 415 9 0.056 4.76 158 0.029 4.35 1020 12 0.159 1.41 0 0.067 1.92 181 15 0.060 3.03 194 0.077 1.80 97 3 0.037 33.33 23 0.036 20.00 37 6 0.068 1.67 32 0.067 1.11 181 2 9 0.141 3.85 0 0.134 2.70 4 12 0.075 1.79 37 0.095 0.90 119 15 0.190 0.01 315 0.707 0.03 119 3 0.011 1.67 207 0.015 9.09 194 6 0.017 5.88 91 0.021 4.17 9 3 9 0.022 6.25 208 0.006 1.67 69 12 0.025 5.88 3 0.017 9.09 49 15 0.024 5.26 18 0.004 9.09 23 3 0.013 7.14 169 0.017 4.76 181 6 0.016 6.67 11 0.013 10.00 14 4 9 0.013 5.56 147 0.003 33.33 111 12 0.047 2.94 119 0.030 3.33 128 15 0.023 1.32 119 0.010 2.70 128 3 0.105 0.83 52 0.088 0.71 60 6 0.005 2.08 0 0.008 0.91 1 5 9 0.002 0.50 23 0.083 0.50 7 12 0.096 0.28 3 0.088 0.42 2 15 0.147 0.28 32 0.054 0.59 28 Table 4.3: Values of the hyperparameter set {-y, c, e}. The (*) denotes the values after WWV was added to the predictor set. The large value of the hyperparameters c(*) reflects the use of nonstandardized predictands in the models. The relative large value of e is due to the large amplitude of the PCs (Figure 2.2). 51 0 C, Chapter 4. Forecast Results (b) Nino3.4 11 70W—i 20W, 5N—5S] (c) Nino3 [i5OW—90W, 5N—5S) 0.9 Ii.i11.1 ead lime (months) (d) Ninol+2 [90W—80W, ON—lOS] 0.9 0.8 0.7 0.6 0.5 L lead time (months) Figure 4.7: Correlation scores before and after WWV was added to the predictor set for the Nino 4, Niño 3.4, Niflo 3, and Niflo 1 + 2 regions. (a) Nino4 1160E—150W, 5N—5S] LR SVR — BNN — 0C, 3 6 9 12 15 52 Chapter 4. Forecast Results C-) 0 0 (c) Nino3 (d) Ninol+2 C-) 0 0 ID Figure 4.8: RMSE scores before and after WWV was added to the predictor set for the Niño 4, Niflo 3.4, Niño 3, and Niño 1 + 2 regions. (b) Nino3.4 C) 0 0 ID 3 6 9 12 15 3 156 9 12 6 9 12 15 lead time (months) lead time (months) 3 53 Chapter 4. Forecast Results Nino3.4 SSTA (1980—2005) (a) 3—month lead : : : : : : — : : Obs. —BNN —LR 1982 1985 1990 1995 2000 2005 (b) 6—month lead : : : : : : : : 1982 1985 1990 1995 2000 2005 (c) 9—month lead : : : : 1982 1985 1990 1995 2000 2005 :(d)12-—mordhlead : : : : : : 1982 1985 1990 1995 2000 2005 (e) 1 5—moith lead : : : 1982 1985 1990 1995 2000 2005 year Figure 4.9: Time series for the Nifio 3.4 region, with BNN (black line), LR (thin line) predictions and observed values (thick grey line). 54 Chapter 4. Forecast Results Nino3.4 SSTA (1980—2005) (a 3—month Idad : : : : . . : : : : — . : : Obs. —SVR —LR 1982 1990 1995 2000 2005 :(by6...onj,)ead: : : : : : : : : I 1982 1985 1990 1995 2000 2005 (c 9—month lead : : : : : : 1982 1985 1990 1995 2000 2005 (d 12.-moith lead . : : : : : : : : 1982 1985 1990 1995 2000 2005 (e 15--mc*ilh lead I I : : : 1982 1985 1990 1995 2000 2005 year Figure 4.10: Time series for the Niño 3.4 region, with SVR (black line), LR (thin line) predictions and observed values (grey line). 55 Chapter 4. Forecast Results 4.11 and 4.12. The advantage in the nonlinear modelling of ENSO shows up at 6 and 9 months lead time in the regions with largest variability (Fig. 4.1 id and f). At these lead times, SVR reveals advantage of nonlinear modelling off the equator along the dateline (Figure 4.12). Adding the WWV index as a predictor effectively reduces the differences in prediction skill between the nonlinear and linear models at long lead times. The prediction advantage at lead times of 9-15 months of BNN is reduced when WWV was added as can be seen by comparing corresponding panels of Figures 4.11 and 4.13. Similarly, the difference in skill between the SVR and linear models is also reduced, however at short lead times of 3 and 6 months the positive differences are enhanced in the eastern and western regions of the domain (Figure 4.14). Compar ing LR against itself after WWV was added, reveals an increase in correlation skill emerging at long lead times of 9-15 months concentrated on the central-eastern equatorial region of the tropical Pacific (Figure 4.15). This increase in skill corre sponds to the region of maximum variability, which seems to indicate that SST in the central-eastern region of the tropical Pacific responds more to changes in the thermocline depth than the western region as formulated in Jin’s oscillator model (1.9). As observed through the forecast of PCi when WWV was the sole predictor (Figure 4.6a) both the linear and nonlinear models indicate a high correlation be tween thermocline depth and SST variations; LR leads the correlation scores at all lead times except for 15 months. 4.2.1 Discussion and Conclusion The forecast results of the 5 leading PCs of SST, indicate evidence of nonlinear structures in the ENSO cycle. This nonlinear structure is mainly present in the sec ond mode of the SST field (Figures 4.1 and 4.5). The central and eastern regions displayed the highest nonlinear tendencies before the WWV is included. The effect of WWV on the dynamics of the SST anomaly along the central-eastern tropical Pacific region appears to be governed by linear dynamics as no advantage was ob served in the forecasts of the nonlinear models when WWV was included (Figures 4.13, 4.14, and 4.15). Stated in [14], the WWV index embodies the large-scale low frequency dynamics. As such, the response of the SST field to the WWV index 56 O T j C &) CD — CD C CD c Z h - _ _ C C - - , p c c CD C. 3 C - CD ç CD - S CD — . — CD CD C - _ I- — CD I’3 — P. 3 P. 3 - . — P. 3 o o rn 0 0 0 0 f l 0 0 ‘ 3I ‘ 3) 0 Z Z CO ‘ 3) 0 Z Z P. 3 — - . P. ) P3 0 0 m 0 0 0 U ) U ) 0 Z Z U ) c J Chapter 4. Forecast Results 1 ON EQ I OS 20S 20N 1 ON EQ lOS (a) SVR (Lead=3 mon) 150E 180 150W 120W 90W (c) SVR (Leod=6 mon) .2u5 -, 150E 180 150W 120W 90W - (e) SVR (Leod=9 mon) iON EQ los 205 20N 1 ON EQ 105 205 2O iON EQ los 150E 180 150W 120W 90W (g) SVR (Leod=12 mon) 150E 180 150W 120W 90W (I) SVR (Leod=15 mon) 15OE IO i5hW 120W 90W (b) SVR — LR 150E 180 150W 120W 90W (d)svR—LR 15’OE 180 150W 120W 90W (f) SVR — LR 20N ....... -_ Figure 4.12: SVR forecast correlation performance and its difference from the linear model for the period 1980-2005 before WWV was included in the predictor set. The contour interval is 0.1 in the right column and 0.05 in the left. iON . . EQ. ‘ Q5 ‘ 1: 105 -. - - E iON 150E 180 150W 120W 90W 3 58 Chapter 4. Forecast Results 20N iON EQ ‘Os. 20S 150E 180 150W 120W 90W 20N (e) BNN(+WWV) (9—mon) ION EQ los 20S T 15OE 180 150W 120W 90W 20N (g) BWW (12—man) iON EQ ‘Os 205 20N iON EQ ‘Os iON EQ los 20S 150E 180 150W 120W 90W 20N - (h) (BNN —LR)(+WWV) Figure 4.13: BNN forecast performance and its difference from the linear model for the period 1980-2005 after WWV was included in the predictor set. The con tour interval is 0.1 in the right column and 0.05 in the left. (a) BNN(+WWV) (3—mon) ION EQ los / . 20N (b) (BNN — LR)(+WWV) 150E 180 150W 120W 90W150E 180 150W 120W 90W (c) BNN(+WWV) (6—mori) p :<— (f) (BNN — LR)(+WWV) 150E 160 150W 120W 96W (I) BNN(+WWV) (15—mon) 150E 180 150W 120W 90W 150E 180 150W 120W SOW 59 Chapter 4. Forecast Results 1 ON EQ los 205 20N 1 ON EQ ws 20S 2Oi 1 ON EQ los (a) SVR(+WWV) (3—mon) 150E 180 150W 120W 90W (c) SVR(+WWV) (6—mon) l5bE 180 150W 120W 90W (e) SVR(-i-WWV) (9—mon) 2O 15OE 180 150W 120W 90W 70N (g)SVR(+WWV) (12—mon) iON EQ 105 205 20N 1 ON EQ los 150E 180 150W 120W 90W (I) SVR(+WWV) (15—rn on) 150E 180 150W 120W 96W (b) (SVR — LR)(+WWV) 150E 180 150W 120W 90W (d) (svR — LR)(+WWV) 150E 180 150W 120W 96W (f) (svR — LR)(+WWV) Figure 4.14: SVR forecast performance and its difference from the linear model for the period 1980-2005 after WWV was included in the predictor set. The con tour interval is 0.1 in the right column and 0.05 in the left. (N. .S. N EQ EQ7 los 150E 180 150W 120W 90W I1 60 00 O -t — - . - - I c < - I — . , I . — (1 C D — o 0 \ Chapter 4. Forecast Results comes with a particular time delay, with a maximum at 6 to 12 months lead time. This observation is in agreement with our results, where the maximum correlation between the SST field and WWV was recorded at 12 months lead time. The forecasts for the 1950-2005 data indicated a marginally nonlinear response in SST at short lead times over certain areas of the domain. These areas are mainly centered off the key Niño regions along the 1800 meridian. This nonlinear struc tures may point to nonlinear effects of the ENSO cycle at higher latitudes [18]. Similar structures appeared in the case of SVR modeling for the later period (1980- 2005) but with higher amplitudes of the correlation scores. In general an increase in the performance of the nonlinear models relative to linear regression was ob served. As argued by [7], changes in the amplitude and spatial distribution of the SST anomalies could signal a changing nature of El Niflo, these changes could be the result of interdecadal variations of the background state, perhaps triggered by global warming. Evidence of these changes could be the recent increased intensity in the warm anomalies, namely the 1982/83 and 1997/98 events. These extreme warm episodes could be responsible for the improved performance of the nonlinear modeling in the shortened time series of 1980-2005 which allowed theses events to stand out. Nevertheless, the increase in correlation skill achieved by isolating the more extreme SST time series anomalies led to only modest advantages in forecasts skill provided by the nonlinear models. Tang et al. (2000) [18] suggested three pos sibilities to explain the limited performance of our nonlinear models. One, the ocean-atmosphere interaction over the tropical Pacific Ocean is essentially linear. Two, the time scales, and data resolution used here are inadequate for extracting the nonlinearities of ENSO as stated in the Cane-Zebiak dynamical model. The third possibility addressed in [18] is that neural network training is not yet optimized, this however is unlikely given the use of SVR which reduces the decisions made during training. The two methods yielded similar results. An important aspect of this study was to investigate the role that warm water volume plays in the development of the SST anomalies. As indicated in [14] the index of WWV embodies the large scale ocean dynamics associated with ENSO. Further confirmation that the the volume of warm water (WWV) above the 20°C isotherm integrated across the basin is linearly related to the ENSO cycle on the 62 Chapter 4. Forecast Results times scales investigated here was obtained [14] [15]. It was also observed that the WWV index alone accounts for as much of the variance as the entire set of predic tors at lead times of 6 months or greater (Figure 4.6a). However, for PC2, WWV alone did not contribute much skill. Likewise, PC3 appears to have no dependence on the WWV index, as correlation skills generally dropped at all lead times for linear and nonlinear regression methods (Figure 4.5). These results suggest that the linear coupled models stated in various versions of the recharge oscillator, e.g. [91, are adequate for modeling the SST and thermocline depth relationship. Fur thermore, the maximum correlation scores between SST and WWV occurred at approximately 12 months lead time which is in agreement with [15]. This increase in correlation scores is concentrated mainly in the eastern region of the tropical Pacific which supports the hypothesis that the dynamical modeling should involve the variables, WWV and TE (SST in the eastern region) with some time lag (Figure 4.15). Forecast of the SST anomalies over the tropical Pacific were conducted using linear and nonlinear methods. Representation of the data in terms of principal components significantly reduced the number of variables needed to be handled, thus avoiding issues associated with large dimensionality. The first 5 PCs of SST account for more than 85% of the total variance of the SST field. Independent forecasts of each of these PCs was carried out using LR, BNN, and SVR. Our results show some general nonlinear features of ENSO, but failed to establish a clear advantage in using BNN or SVR over LR in the forecast of SST, particularly in the analysis of the 1950-2005 period. The nonlinear features showing up in the forecasts of SST were of second order in the variance, namely PC2, these are located along the regions of largest variability (central-eastern tropical Pacific). This nonlinear advantage was reduced when WWV was added to the predictor set. We concluded that ENSO, at interseasonal times scales is essentially linear, the nonlinear features account for only a small percentage of the variance. As pointed out by [18], other nonlinear elements of ENSO could be contained in the seasonal cycle. However, the seasonal cycle was eliminated in this analysis. 63 Chapter 4. Forecast Results 4.3 Further Research The forecast model development presented here could be extended to a third non linear model known as Gaussian Processes. Gaussian Processes is another kernel method, which extends the use of kernels to the realm of probabilistic forecast in the Bayesian setting. There are many aspects of ENSO left uninvestigated. There is strong evidence that many processes associated with this phenomenon exhibit nonlinear behavior, for example the locking of El Niño events with the seasonal cycle are difficult to explain with linear methods, and the asymmetries between El NiñofLa Nina am plitudes. Low order chaotic forcing has been proposed as a mechanism for locking ENSO to the seasonal cycle [19]; asymmetries are attributed to nonlinear dynami cal heating [2]. A statistical nonlinear analysis could be applied to investigate the degree, if any, these processes have on the ENSO cycle. Furthermore, Tang et al. [181 proposed that nonlinear methods should also be used to study the response of the extratropics to the tropical SST anomalies. Of particular interest is to study the effect of high frequency atmospheric forc ing on ENSO variability as an extension of the present project. McPhaden et al. (2006) used a simple linear regression model with an index of the Madden Julian Oscillation (MJO) and the (WWV) index as predictors for the tropical Pacific SST. This simple two predictor model reportedly accounts for -‘ 60% of the total SST variance at 2-3 season lead times [14]. A natural extension of this work would be to apply the nonlinear methods, BNN and SVR to this system to investigate if high frequency forcing characterized by the MJO index has a nonlinear influence on low frequency ENSO dynamics. 64 Bibliography [1 J Allen, R.; Lindsey, J. and Parker, D. (1996). El Niflo Southern Oscillation and climatic variability. CSIRO Plublishing. [2] An, S.-I. and 3m F.-F. (2004). Nonlinearity and Asymmetry of ENSO. Journal of Climate, 17(12):2399—2412. [3] Bishop, C.M. (2006). Pattern Recognition and Machine Learning. Library of Congress. [4] Bjerknes, J. (1969). Atmospheric Teleconnection from the Equatorial Pacific. Monthly Weather Review, 97(3): 163—172. [5] Chang, C.-C. and Lin, C.-J. (2001). LIBSVM: a library for support vec tor machines. Software available at http: / /www. csie . ntu. edu. tw/ ‘cj lin/libsvm. [6] Clarke, A.J.; Van Gorder, S.; Colantuouno, G. (2007). Wind stress curl and ENSO discharge/recharge in the equatorial Pacific. Journal of Physical Oceanography, 37(4): 1077—1091. [7] Fedorov, A.V. and Philander, S.G. (2000). Is El Niño Changing? Science, 288(3): 1997—2002. [8] Hsieh, W.W. (2009). Machine Learning Methods in Environmental Sciences — Neural Networks and Kernels. Cambridge University Press. [9] Jin, F.-F. (1997). An equatorial ocean recharge paradigm for ENSO. Part I and II: Conceptual model. Journal of the Atmospheric Sciences, 54(7):811—847. [10] Kalnay, M.E.; et at. (1996). The NCEP/NCAR reanalysis project. Bulletin of the American Meteorological Society, 77(3):437—47 1. 65 Bibliography [11] Katz, Richard W. (2002). Sir Gilbert Walker and a connection between El Niflo and statistics. Statistical Science, 17(1):97—102. [12] Latif, M.; Bamett, T.P.; Cane, M.A.; Flugel, M.; Graham, N.E.; von Storch, H.; Xu, 3.-S. and Zebiak, S.E. (1994). A review of ENSO prediction studies. Earth and Environmental Science, 9(4—5):167—179. [13] Marzban, C. (2004). My Neural Network short course presented at the American Meteorological Society Conference. http: / /www. nhn. ou. edu/ marzban. [14] McPhaden, M.J.; Zhang, X.; Hendon, H.H. and Wheeler. M.C. (2006). Large scale dynamics and MJO forcing of ENSO variability. Geophysical Review Letters, 33:L 16702. [15] Meinen, C.S. and McPhaden, M.J. (2000). Observation of warm water vol ume changes in the equatorial Pacific and their relationship to El Niño and La Nina. Journal of Climate, 13(20): 3551—3559. [16] Philander, S .G. (1990). El Niño, La Nina, and the Southern Oscillation, vol ume 46 of International Geophysics Series. Academic Press, mc, San Diego; California, first edition. [17] Smith, T.M. and Reynolds, R.W. (2004). Improved extended reconstruction of SST. Journal of Climate, 17(12):2466—2477. [181 Tang, B.; Hsieh, W.W.; Monahan, A.H. and Tangang, FT. (2000). Skill Com parisons between Neural Networks and Canonical Correlation Analysis in Pre dicting the Equatorial Pacific Sea Surface Temperatures. Journal of Climate, 13(1):287—293. [19] Tziperman, E.; Cane, M.A. and Zebiak, S.E. (1995). Irregularity and locking to the seasonal cycle in an ENSO Prediction model as explained by the quasi periodicity route to chaos. Journal of the Atmospheric Sciences, 52(3):293—306. [20] van den Dool, A.G. et.al (1994). Long-lead seasonal forecast: where do we stand? Bulletin of the American Meteorological Society, 75(1 1):2097—21 14. 66 Bibliography [211 von Storch, H. and Zwiers, EW. (1999). Statistical Analysis in Climate Re search. Cambridge University Press. [22] Wang, W. and McPhaden, M.J. (2000). The surface-layer heat balance in the equatorial Pacific ocean. part II: interannual variability. American Meteorologi cal Society, 30(1 1):2989—3008. [23] Wu, A.; Hsieh, W.W. and Tang, B. (2006). Neural network forecasts of the tropical Pacific sea surface temperatures. Neural Networks, 19(2):145—154. [24] Wyrtki, K. (1975). El Niño—the dynamic response of the equatorial Pacific ocean to atmospheric forcing. Journal of Physical Oceanography, 5(4):572— 584. 67 Appendix A Time Series for El Niño 4,3.4,3, and 1+2 68 Appendix A. Time Series for El Niño 4, 3.4, 3, and 1÷2 Nino4 SSTA (1980—2005) (a> 3—month lead 1982 1985 1990 1995 2000 2005 2 (b) 6—month lead 1982 1985 1990 1995 2000 2005 2 (c 9—month lead : 1982 1985 1990 1995 2000 2005 2 (d) 12—month lead : 1982 1985 1990 1995 2000 2005 2 (e)15—monthlead 1982 1985 1990 1995 2000 2005 year Figure A. 1: Time series for the Niño 4 region, showing the BNN (black line), and LR (thin line) predictions and the observed values (thick grey line). 69 Appendix A. Time Series for El Niflo 4, 3.4, 3, and 1+2 Nino4 SSTA (1 980—2005) (+WWV) 2 :(a)34flOflthIead Obs. BNN LR 1982 1985 1990 1995 2000 2005 (b6_rnonthIéad: : 1982 1985 1990 1995 2000 2005 2 (c 9—month lead : 1982 1985 1990 1995 2000 2005 2 (d)12_mQ,thIead:::: 1982 1985 1990 1995 2000 2005 2 : (e15.—monthIead 1982 1985 1990 1995 2000 2005 year Figure A.2: Time series for the Niflo 4 region, showing the BNN and LR predic tions. WWV is included in the predictor set. 70 Appendix A. Time Series for El Niflo 4, 3.4, 3, and 1÷2 Nino4 SSTA (1980—2005) 2 (a) 3—month lead ; : 1982 1985 1990 1995 2000 2005 (b) 6—month lead : : 1982 1985 1990 1995 2000 2005 (c) 9—month lead 1982 1985 1990 1995 2000 2005 (d) 12—month lead : : : : : : 1982 1985 1990 1995 2000 2005 (e) 15—month lead 1982 1985 1990 1995 2000 2005 year Figure A.3: Time series for the Niño 4 region, showing the SVR and LR predic tions. 71 Appendix A. Time Series for El Niño 4, 3.4, 3, and 1+2 Nino4 SSTA (1980—2005) (+WWV) (a) 3—month lead : 1982 1985 1990 1995 2000 2005 (b) 6—month lead 1982 1985 1990 1995 2000 2005 2 (c 9—month lead : 1982 1985 1990 1995 2000 2005 2 (d) 12—month lead : : : 1982 1985 1990 1995 2000 2005 (e) 1 5—month lead : : : : : : : : 1982 1985 1990 1995 2000 2005 year Figure A.4: Time series for the Niflo 4 region, showing the SVR and LR predic tions. WWV is included in the predictor set. 72 Appendix A. Time Series for El Niño 4, 3.4, 3, and 1+2 Nino3.4 SSTA (1980—2005) (÷WWV) (a) 3—month lead : : : : :.: : — : : : Obs. BNN —La 1982 1985 1990 1995 2000 2005 (b) 6—month lead : : : : 1982 1985 1990 1995 2000 2005 :(cy9_monthlead: : : : : : : I : 1982 1985 1990 1995 2000 2005 (d 1 2—month lead : : : S : : 1982 1985 1990 1995 2000 2005 (e) 15—month lead S : 1982 1985 1990 1995 2000 2005 year Figure A.5: Time series for the Nino 3.4 region, showing the BNN and LR pre dictions. WWV is included in the predictor set. 73 Appendix A. Time Series for El Niño 4, 3.4, 3, and 1+2 Nino3.4 SSTA (1980—2005) (+WWV) (a)3_monhlèad : : : : : : —2 Obs —SVR —LR 1982 1985 1990 1995 2000 2005 (b 6—monTh lead : : : : : : : : 1982 1985 1990 1995 2000 2005 .(cy9—nlonthlead : : : : : : 1982 1985 1990 1995 2000 2005 (d)12——rnonthIead : : : : : : : : 1982 1985 1990 1995 2000 2005 :(e)15mci,thIead : : : : : : : : 1982 1985 1990 1995 2000 2005 year Figure A.6: Time series for the Niño 3.4 region, showing the SVR and LR predic tions. WWV is included in the predictor set. 74 Appendix A. Time Series for El Niño 4, 3.4, 3, and 1+2 Nino3 SSTA (1 980—2005) —2 Obs — BNN LB 1982 1985 1990 1995 2000 2005 month lead 1982 1985 1990 1995 2000 2005 2month lead 1982 1985 1990 1995 2000 2005 1982 1985 1990 1995 2000 2005 çmthIead 1982 1985 1990 1995 2000 2005 year Figure A.7: Time series for the Niño 3 region, showing the BNN and LR predic tions. 75 Appendix A. Time Series for El Niño 4, 3.4, 3, and 1+2 Nino3 SSTA (1980—2005) (+WWV) .. —2 ois — BNN LR 1982 1985 1990 1995 2000 2005 month lead 1982 1985 1990 1995 2000 2005 1982 1985 1990 1995 2000 2005 . 1982 1985 1990 1995 2000 2005 çonthlead 1982 1985 1990 1995 2000 2005 year Figure A.8: Time series for the Niño 3 region, showing the BNN and LR predic tions. WWV is included in the predictor set. 76 Appendix A. Time Series for El Niño 4, 3.4, 3, and 1+2 Nino3 SSTA (1980—2005) —2 . Obs. — SVR LR 1982 1985 1990 1995 2000 2005 2onth lead 1982 1985 1990 1995 2000 2005 2m0nhh lead 1982 1985 1990 1995 2000 2005 2mth1ea 1982 1985 1990 1995 2000 2005 çrnonthlead 1982 1985 1990 1995 2000 2005 year Figure A.9: Time series for the Niño 3 region, showing the SVR and LR predic tions. 77 Appendix A. Time Series for El Niño 4, 3.4, 3, and 1+2 Nino3 SSTA (1980—2005) (+WWV) 2pçheJ —2 Obs —SVF —LR 1982 1985 1990 1995 2000 2005 2ad 1982 1985 1990 1995 2000 2005 1982 1985 1990 1995 2000 2005 2hlead 1982 1985 1990 1995 2000 2005 2çthead 1982 1985 1990 1995 2000 2005 year Figure A. 10: Time series for the Niño 3 region, showing the SVR and LR predic tions. WWV is included in the predictor set. 78 Appendix A. Time Series for El Niño 4, 3.4, 3, and 1÷2 Ninol+2 SSTA (1980—2005) —3 Obs. —BNN —LR 1982 1985 1990 1995 2000 2005 4 .:()_çQfl lead; ; ; V 1982 1985 1990 1995 2000 2005 1982 1985 1990 1995 2000 2005 4çrnthie 1982 1985 1990 1995 2000 2005 1982 1985 1990 1995 2000 2005 year Figure A. 11: Time series for the Niflo 1 + 2 region, showing the BNN and LR predictions. 79 Appendix A. Time Series for El Niño 4, 3.4, 3, and 1+2 Ninol+2 SSTA (1980—2005) (+WWV) — Obs. — BNN LR 1982 1985 1990 1995 2000 2005 1982 1985 1990 1995 2000 2005 1982 1985 1990 1995 2000 2005 1982 1985 1990 1995 2000 2005 1982 1985 1990 1995 2000 2005 year Figure A. 12: Time series for the Niño 1 + 2 region, showing the BNN and LR predictions. WWV is included in the predictor set. 80 Appendix A. Time Series for El Nino 4, 3.4, 3, and 1+2 Ninol+2 SSTA (1980—2005) 3 Obs —SVR —L 1982 1985 1990 1995 2000 2005 1982 1985 1990 1995 2000 2005 1982 1985 1990 1995 2000 2005 1982 1985 1990 1995 2000 2005 1982 1985 1990 1995 2000 2005 year Figure A. 13: Time series for the Niño 1 + 2 region, showing the SVR and LR predictions. 81 Appendix A. Time Series for El Niiio 4, 3.4, 3, and 1+2 Ninol+2 SSTA (1980—2005) (÷WWV) Obs —SVR —LR 1982 1985 1990 1995 2000 2005 1982 1985 1990 1995 2000 2005 1982 1985 1990 1995 2000 2005 1982 1985 1990 1995 2000 2005 1982 1985 1990 1995 2000 2005 year Figure A.14: Time series for the Niño 1 + 2 region, showing the SVR and LR predictions. WWV is included in the predictor set. 82
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Forecasts of tropical Pacific sea surface temperatures...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Forecasts of tropical Pacific sea surface temperatures by neural networks and support vector regression Aguilar-Martinez, Silvestre 2008
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Forecasts of tropical Pacific sea surface temperatures by neural networks and support vector regression |
Creator |
Aguilar-Martinez, Silvestre |
Publisher | University of British Columbia |
Date Issued | 2008 |
Description | Nonlinear and linear regression models were developed to forecast the sea surface temperature anomalies (SSTA) across the tropical Pacific ocean. The methods used were, Bayesian neural networks (BNN), support vector machines for regression (SVR) and linear regression (LR). The predictors of the models were a particular combination of the principal components of sea level pressure (SLP) and sea surface temperature (SST) data at lead times ranging from 3 to 15 months. Two data sets corresponding to the time periods, 1950-2005 and 1980-2005 were independently studied. The later data set contained the standardized volume of warm water (WWV) integrated across the tropical Pacific basin as an extra predictor. The seasonal cycle was removed from both data sets prior to their introduction into the regression models. Various graphical and statistical tools were used to compare the performance of these models, in particular the correlation scores. The results reveal that although there appeared to be nonlinear elements of the SSTA response to the previous SLP and SST conditions, this response is of second order. The large amplitude variance represented by the first principal component (PC) of the SST field can be adequately expressed as a linear response. Consistently in both data sets, the advantages in the nonlinear modelling manifested mainly through the second PC of the SST variability. The spacial distribution of correlation skills between the predicted and observed SST fields corresponding to the BNN and SVR models show some nonlinear structures in various parts of the domain. Despite some differences in the temporal and spacial distribution of these correlation skills, there was no decisive advantage between using BNN or SVR. Adding the WWV to the predictor set reduced the marginal advantage of the nonlinear methods over the eastern region of the tropical Pacific, in particular at lead times of 9-15 months. These results generally agree with linear recharge oscillator models, and the hypothesis of a linear relationship between thermocline depth and SST in the eastern region. |
Extent | 2851110 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-04-27 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0067175 |
URI | http://hdl.handle.net/2429/7577 |
Degree |
Master of Science - MSc |
Program |
Physics |
Affiliation |
Science, Faculty of Physics and Astronomy, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2009-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2009_spring_aguilar-martinez_silvestre.pdf [ 2.72MB ]
- Metadata
- JSON: 24-1.0067175.json
- JSON-LD: 24-1.0067175-ld.json
- RDF/XML (Pretty): 24-1.0067175-rdf.xml
- RDF/JSON: 24-1.0067175-rdf.json
- Turtle: 24-1.0067175-turtle.txt
- N-Triples: 24-1.0067175-rdf-ntriples.txt
- Original Record: 24-1.0067175-source.json
- Full Text
- 24-1.0067175-fulltext.txt
- Citation
- 24-1.0067175.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}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0067175/manifest