Analysing time-varying trends in stratospheric ozone time series using the state space approach

We describe a hierarchical statistical state space model for ozone profile time series. The time series are from satellite measurements by the Stratospheric Aerosol and Gas Experiment (SAGE) II and the Global Ozone Monitoring by Occultation of Stars (GOMOS) instruments spanning the years 1984–2011. Vertical ozone profiles were linearly interpolated on an altitude grid with 1 km resolution covering 20– 60 km. Monthly averages were calculated for each altitude level and 10 wide latitude bins between 60 ◦ S and 60 N. In the analysis, mean densities are studied separately for the 25–35, 35–45, and 45–55 km layers. Model variables include the ozone mean level, local trend, seasonal oscillations, and proxy variables for solar activity, the Quasi-Biennial Oscillation (QBO), and the El Niño–Southern Oscillation (ENSO). This is a companion paper to Kyrölä et al.(2013), where a piecewise linear model was used together with the same proxies as in this work (excluding ENSO). The piecewise linear trend was allowed to change at the beginning of 1997 in all latitudes and altitudes. In the modelling of the present paper such an assumption is not needed as the linear trend is allowed to change continuously at each time step. This freedom is also allowed for the seasonal oscillations whereas other regression coefficients are taken independent of time. According to our analyses, the slowly varying ozone background shows roughly three general development patterns. A continuous decay for the whole period 1984–2011 is evident in the southernmost latitude belt 50–60 ◦ S in all altitude regions and in 50–60 ◦ N in the lowest altitude region 25–35 km. A second pattern, where a recovery after an initial decay is followed by a further decay, is found at northern latitudes from the equator to 50 ◦ N in the lowest altitude region (25–35 km) and between 40 ◦ N and 60 N in the 35–45 km altitude region. Further ozone loss occurred after 2007 in these regions. Everywhere else a decay is followed by a recovery. This pattern is shown at all altitudes and latitudes in the Southern Hemisphere (10–50 ◦ S) and in the 45– 55 km layer in the Northern Hemisphere (from the equator to 40 N). In the 45–55 km range the trend, measured as an average change in 10 years, has mostly turned from negative to positive before the year 2000. In those regions where the “V” type of change of the trend is appropriate, the turning point is around the years 1997–2001. To compare results for the trend changes with the companion paper, we studied the difference in trends between the years from 1984 to 1997 and from 1997 to 2011. Overall, the two methods produce very similar ozone recovery patterns with the maximum trend change of 10 % in 35–45 km. The state space method (used in this paper) shows a somewhat faster recovery than the piecewise linear model. For the percent change of the ozone density per decade the difference between the results is below three percentage units.


Introduction
Time series constructed from satellite remote sensing observations provide important information about variability and trends in atmospheric chemical composition.Many satellite time series provide global coverage of the measurement and some of the data sets run from the 1980s.The analysis of trends, both natural and anthropogenic, is complicated by natural variability and forcing affecting stratospheric chemical compositions.In this study, the recovery of stratospheric ozone from the depletion caused by Published by Copernicus Publications on behalf of the European Geosciences Union.chlorofluorocarbon (CFC) compounds is studied using a statistical time series model.
Slow background changes in stratospheric ozone are easily masked by both seasonal and irregular natural variabilities.Thus, the requirements are stringent for the stability of ozone observations.Self-calibrating occultation instruments are good candidates for such a task.The observations analysed in this work consist of satellite measurements by the SAGE II and GOMOS instruments, operational during 1984-2005and 2002-2012, respectively.The data set used here spans the years 1984-2011.Vertical ozone profiles were linearly interpolated on an altitude grid with 1 km resolution covering 20-60 km.Monthly averages were calculated for each altitude level and 10 • wide latitude bins between 60 • S and 60 • N. Combining the observations from different instruments having different measurement principles is a challenge.Kyrölä et al. (2013) explain the data set and its construction in more detail.Here the analysis is done with both the original 1 km vertical spacing and by calculating mean densities over 10 km intervals.
There is a wealth of literature concerning the analysis of atmospheric time series.A good reference to stratospheric ozone time series regression analysis is SPARC (1998).A recent study that reviews the challenges and problems in trend analysis of climatic time series was published by Bates et al. (2012), and a general trend analysis reference is Chandler and Scott (2011).For state space and functional analysis of atmospheric time series of similar type to that performed here, see Lee and Berger (2003) and Meiring (2007).
This paper studies the feasibility and practical implementation of a state space approach for atmospheric time series analysis by defining a dynamic linear model (DLM) for stratospheric ozone time series."Dynamic" means here that the regression coefficients can evolve in time.This makes it possible to describe and analyse smooth changes in the average background behaviour of ozone.Model variables include the ozone mean level, local trend, seasonal oscillations, and proxy variables for solar activity, the Quasi-Biennial Oscillation (QBO), and the El Niño-Southern Oscillation (ENSO).We do not claim novelty in the presented methods themselves, but argue that they should be more extensively applied in the analysis of climatic time series and provide a simple framework for time series analyses that can be generalized to more comprehensive studies.In this paper, we describe the necessary steps for applying the methods.
A typical feature in atmospheric time series is that they are not stationary but exhibit both slowly varying and abrupt changes in the distributional properties.These are caused by external forcing such as changes in the solar activity or volcanic eruptions.Further, the data sampling is often nonuniform -there are data gaps, and the uncertainty of the observations can vary.When observations are combined from various sources there will be instrument and retrieval method related biases.The differences in sampling lead to uncertainties, too.Straightforward linear regression analysis leaves the model residuals correlated as not all variability can be explained by a static linear structure.Usually this is compensated by allowing some correlation structure to the model observation error by using, e.g. an autoregressive model.If the residual correlation is not accounted for, the model uncertainty analyses are misleading.A simple autoregressive process can explain some of the unmodelled systematic variations by correlated noise, again confusing the analyses.In conclusion, much care in interpretation is needed for those standard classical statistical time series methods that require stationarity, such as the ARIMA (autoregressive integrated moving average) approach.The more general approach discussed in this paper makes use of dynamic linear models and Kalman filter type sequential estimation algorithms.
State space models, sometimes called hidden Markov models or structural time series models, are well known and documented in time series literature, e.g.Harvey (1990); Hamilton (1994); Migon et al. (2005).Modern computationally oriented references are Durbin and Koopman (2012) and Petris et al. (2009).Here, we review the basic properties relevant to the analysis of atmospheric ozone time series data and explain the necessary steps to fit the model to monthly time series observations and how to assess the uncertainties in the trend estimation.
The structure of this paper is the following.The data sets and the statistical model are described in Sect. 2. Results of the statistical time series analyses are given in Sect. 3 and the paper ends with discussion and conclusions in Sect. 4. Appendix A contains mathematical and computational details.

Ozone time series from satellite observations
We use a combination of two ozone data sets.The first consists of solar occultation measurements of ozone in the stratosphere and lower mesosphere from the SAGE II instrument (Chu et al., 1989(Chu et al., ) operational during 1984(Chu et al., -2005.The second is the GOMOS instrument (Bertaux et al., 2010) that measured ozone in the stratosphere, mesosphere and lower thermosphere during 2002-2012 using stellar occultations.The individual data sets have been homogenized to form a combined time series from 1984 to 2011.The stability of the SAGE II and GOMOS instruments, the construction of the combined time series, data screening, bias correction, and other issues are discussed in more detail in the companion paper by Kyrölä et al. (2013).Four proxy time series are used for the solar flux, Quasi-Biennial Oscillations and ENSO.Except for ENSO, these proxies are the same as used in the companion paper.For ENSO we use the Multivariate ENSO Index (MEI) from NOAA1 .

Statistical time series model
A general linear state space model with Gaussian errors can be written with an observation equation and a state evolution equation as where y t is a vector of length p containing the observations and x t is a vector of length q of unobserved states of the system at time t.The state variables are used to describe the various components of the time series model, such as mean level, trend, seasonality and the effect of proxy variables.Matrix F t is the observation operator mapping the unobserved states to the observations and matrix G t is the model evolution operator giving the dynamics of the states.In this basic formulation the uncertainties v t and w t are assumed to be Gaussian, with observation uncertainty covariance V t and model error covariance W t .Above, N p (0, V t ) stands for p-dimensional Gaussian distributions, with vector of zeros as mean and V t as the p×p covariance matrix.The time index t will go from 1 to N, the length of the time series to be analysed.In the following, the matrices defining the model will mostly be time invariant, i.e.G t = G, etc., and we will usually drop the time subscript, still retaining it in general formulas that are not specific to this particular time series application.
In this work, we use a DLM to explain variability in the ozone time series with four components: smooth locally linear trend, seasonal effect, effect of forcing via proxy variables, and noise that is allowed to have autoregressive correlation.All components are built using the state space approach.
To describe the trend we start with a simple local level and locally linear trend model that has two hidden states x t = [µ t , α t ] T , where µ t is the mean level and α t is the change in the level from time t to time t + 1.In addition we need stochastic terms for the observational error and for the allowed change in the dynamics of the trend and the level.These are defined by Gaussian " " terms below.The system can be written by equations In terms of the state space equations (1) and (2) this becomes Depending on the choice of variances σ 2 level and σ 2 trend , this will define a smoothly varying background level of the data series that is used to infer changes in atmospheric ozone.
Most atmospheric series exhibit seasonal variability.The seasonality can be modelled using harmonic functions.If the number of cyclic components is s, the full seasonal model has s/2 harmonics.For the kth harmony, with k = 1, . . ., s/2, we need to add two state variables.With monthly data we have s = 12 and the corresponding blocks of the model and observation matrices are Here the state equation matrices are independent of time index t and we have used subscript k to stand for the harmonic component.
where G seas(k) is defined in Eq. ( 7) and does not depend on time t (e.g.Petris et al., 2009, Sect. 3.2.3).In our case the seasonality can be adequately explained by two harmonics, e.g. by yearly and half-a-year variation (k = 1 and k = 2), which will increase the number of hidden states to be estimated by four.In addition we need to define the error covariance matrix W seas for the allowed time-wise variability in the seasonal components.We note that there is no need for a separate seasonal observation error matrix V seas , as V t is used only in the observation equation ( 1), not in the model equation ( 2) where the seasonal cycle is defined.
In the previous definitions, the state operator G, the observation operator F and the model error covariance W have been time invariant.The observation uncertainty covariance V t defined in Eq. ( 6) is, in our case, time dependent and it will contain the known observation uncertainties.The inclusion of auxiliary proxy variables is done by augmenting the observation matrix F t , making it time dependent.In the following, the stratospheric ozone analysis will utilize four proxy time series explaining parts of the natural variability: one for the solar flux, two proxy variables for QBO, and one for ENSO.This is achieved by adding the following components into the system matrices: where z 1,t , z 2,t , z 3,t and z 4,t contain the values of the four proxy series at time t.With positive error variances σ 2 zi this is an extension of linear regression analysis into one with timevarying coefficients.In our analysis, for simplicity, we set all the elements of the model error covariance matrix W proxy to zero to obtain time-invariant regression coefficients for the proxy variables.
To allow autocorrelation in the residuals we use a firstorder autoregressive model (AR(1)).This is similar to the Cochrane-Orcutt correction in classical multiple regression (see e.g.Hamilton, 1994) used by Kyrölä et al. (2013).However, in the DLM approach we can estimate the autocorrelation coefficient and the extra variance term together with the other model parameters, not by a separate iteration, as needed, e.g. by the Cochrane-Orcutt method.A first-order autoregressive model for a state component η t can be written as η t = ρη t−1 + AR , with AR ∼ N(0, σ 2 AR ), where ρ is the AR coefficient and σ 2 AR is usually called the innovation variance in classical time series literature.In DLM form, we simply have and both ρ and σ 2 AR can be estimated from the observations.The next step in DLM model construction is the combination of the selected individual model components into larger model evolution and observation equations.For the model evolution matrix G and the observation operator F t we have In our analyses, the state vector x t has a total of 11 components: where µ t is the local mean, α t is the local trend, i.e. the local change of the mean, u t,k and u * t,k are the states representing the two seasonal harmonics, the β are the regression coefficients for the four proxy variables, and η t is an extra state for the autoregressive component that "remembers" the value of the observation from the previous time step.In this study, the regression coefficients (β) are assumed not to depend on time t.The full model error covariance matrix W is a diagonal matrix with the corresponding variances at the diagonal, where the local level variance and the proxy regression coefficient variances have been set to zero and the four seasonal variances set to be all equal, to correspond to the simplification assumed in the analyses.Lastly, as already stated, the observation error covariance matrix V t , will be a 1×1 matrix and equals the known observation uncertainty: V t = σ 2 obs(t) .Next, the analysis will proceed to the specification of the variance parameters and other parameters in the model formulation (e.g. the AR coefficient ρ above), and to the estimation of the model states by state space methods.

Model parameter estimation
We have two kinds of unknowns, the model state variables x t , one vector for each time t, and the auxiliary parameters that define the model error covariance matrix W and the system evolution matrix G.At first sight, this might seem to be a vastly underdetermined system, as we have several model parameters for each single observation.However, by the sequential nature of the equations, we can estimate the states by standard recursive Kalman filter formulas.
Implicitly assumed in the state space equations (1) and ( 2) is that the state at time t is statistically conditionally independent of the history given the previous state at time t − 1.When the model equation matrices are known, this Markov property allows sequential estimation of the states given the observations by famous Kalman formulas (see e.g.Rodgers, 2000).We can use a Kalman filter for one-stepahead prediction of the state and a Kalman smoother for the marginal probability distribution of the state at time t given the whole time series of all observations (t = 1, . . ., N ).Marginal means that the uncertainty of the states at all other times than t has been accounted for.In time series applications, the Kalman filter output can be used to calculate the model likelihood function, needed in the estimation of the auxiliary parameters, and the Kalman smoother provides an efficient algorithm to estimate the model states and to decompose time series into parts given by the model formulation.Thus, this high-dimensional problem is computationally not much more intensive than a classical static multiple linear regression analysis.Furthermore, in the linear Gaussian case, the probability distributions provided by the Kalman formulas are exact, not approximations.Non-linear models can be approached by linearization of the state equations, and non-Gaussian error models by, e.g.particle filter algorithms (Doucet et al., 2001).More details on the DLM computations can be found in Appendix A.
Next, we consider the model error covariance matrix W. If we set all model error variances to zero it will change the DLM model into an ordinary, non-dynamic, multiple linear regression model.By using non-zero variances we can fit a smoothly varying mean level, and the smoothness can be controlled by the size of the variances.A simplification done here is that we assume W to be diagonal, and even some of the diagonal elements are set to zero.In our case non-zero elements are the variability of the trend, σ2 trend , and the variances of seasonal variation, σ 2 seas(k) , i.e. we will set the variance of the level and the four proxy variables to zero.Our motivation for excluding dynamic variation from the proxies is to use an as simple as possible dynamic model that would have similar properties as the piecewise linear model in the companion paper.By studying dynamic changes in the trend and seasonal variation we will then either validate the static linear approach or show that it is not appropriate.
A common procedure to estimate the elements of the model error matrix W is based on the maximum likelihood method using the likelihood function provided by the Kalman filter.After the estimation, the obtained values could be plugged into the system equations as known constants.However, this plug-in method neglects the uncertainty in the estimates.Instead, we will use an alternative method based on Bayesian analysis (Gamerman, 2006;Petris et al., 2009;Särkkä, 2013) and outline it shortly below in Sect.2.4.

Markov chain Monte Carlo analysis for model parameters
The Markov chain Monte Carlo (MCMC) method provides an algorithm to draw samples from the Bayesian posterior uncertainty distribution of model parameters given the likelihood function and the prior distribution (Gamerman, 2006).
As the Kalman filter likelihood provides the likelihood of the auxiliary model parameters, marginalized over the unknown model states, we can use it to draw samples from the marginal posterior distribution of these parameters.We will use an adaptive Metropolis algorithm of Haario et al. (2006) for the three unknown variance parameters (σ trend , σ seas , σ AR ) in the matrix W in Eq. ( 13) as well as for the autoregressive parameter ρ in the system evolution matrix G (Eq. 10).After obtaining this sample from the posterior distribution, we use sampled parameter values, one by one, to simulate realizations of the model states x 1:N using the Kalman simulation smoother (e.g.Petris et al., 2009).This allows us to account for both the parameter and the model state uncertainties in the trend analysis.Again, more computational details are given in Appendix A.
For Bayesian analysis of the unknown model parameters, prior distributions for the parameters must be specified.We want the variances in the matrix W to reflect our prior knowl-edge on the assumed variability in the processes captured by the observations.As noted by Gamerman (2006, Sect. 2.5.), dynamic linear models offer intuitive means of providing qualitative prior information in the form of the model equations and quantitative information by prior distributions on variance parameters.By the estimation procedure we aim at finding variance parameters that are consistent with the given observation uncertainty, i.e. the model can predict the observations within their accuracy.This means that the scaled prediction residuals that are used defining the likelihood function should behave like an independent Gaussian random variable.We can assess these assumptions by different residual analysis diagnostics.
As we are effectively looking for slowly varying trends in the data, we will set prior constraints to variance parameters to reflect this.For example, we might assume that the change within a month in the background level is on the average some percentage of the overall time series mean.The estimation procedure will then divide the observed variability into model components (level, trend, seasonality) in proportions that reflect the prior choices.The standard model diagnostic tools, such as autocorrelation analysis and normal probability plots, can be used to reveal possible discrepancies in the model assumptions that have to be considered.As the model residuals are calculated from one-step predictions, the diagnostics will reveal both over-fit and a lack of fit.
The model error covariance W (Eq. 13) is time invariant, with nonzero diagonal terms for the trend parameter, a common value for the four parameters defining the variability in the seasonal components, and one value for the variance defining the autoregressive component.For MCMC estimation of the variance parameters, we use log-normal prior distributions for the corresponding standard deviations.The motivation for this is that the standard deviations are positive by definition so logarithmic scale is natural and commonly used and it allows prior specification in meaningful units.If a random variable x follows a log-normal distribution logN (µ, σ 2 ), then log(x) follows the standard Gaussian distribution N (µ, σ 2 ).Because of this transformation, it is more intuitive to work with a log mean parameter, as in logN (log(µ), σ 2 ), where µ is now the approximate mean of the untransformed x and the scale parameter of the lognormal distribution σ can be interpreted as an approximation of the relative standard deviation of the original variable in question.The prior distribution for the standard deviation of the monthly level change σ trend was set to be log-normal with the log mean equal to 1/12 % of the average level of the ozone observations and having scale parameter 1.This scale corresponds to a relatively wide (130 % SD) 2 prior uncertainty for σ trend .All the σ parameters describe allowed change in one time interval of the model and data (month), i.e. the unit for σ trend is 1 cm −3 .For the seasonal standard deviation, σ seas , we set the log mean parameter to 1 % of the overall variability in the observed ozone values and the scale parameter to 2 (corresponding to 730 % SD).For the standard deviation of the AR(1) component (Eq.10), we use 30 % of the observed variability for the log mean and similarly 2 for the scale.For the autoregressive coefficient ρ the prior was Gaussian N(0.45, 0.5) truncated to be in the interval [0, 1].Table 1 collects the prior specifications.The results for the trend analyses or the estimated parameter values themselves are not sensitive to the chosen priors.The MCMC for each separate fit was run using 10 000 simulation steps, and the resulted parameter samples (chains in MCMC terminology) were checked for convergence of the algorithm.The histograms of the marginal posterior samples of the parameters with the prior probability density functions overlaid are shown in Fig. 1 for one zonal band and altitude region.The prior for the trend variability would allow the trend to change    2013) is drawn with a solid green line.The second panel shows model residuals, scaled by the estimated residual standard deviation, i.e. the variability in the y axis direction should be approximately standard Gaussian N (0, 1).The other panels show means and 95 % probability envelopes (grey shading) of the DLM components related to seasonality and to the proxy variables -solar, QBO and ENSO -i.e. the original proxy variables multiplied by the estimated regression coefficients.For the proxy variables, the y axis scaling is the percentage of the total ozone variability, i.e. the ozone observations were scaled by their standard deviation, individually for each data set analysed.At 20 km we see a large number of missing observations, some possible outliers, and a strong effect of the ENSO index.There is only a slight difference between DLM and linear models, however, although the latter does not include ENSO.
by about same amount as it does in linear trend analyses -or let it not change if the data do not support it -but, also, should not capture too fine changes due to natural variability that could be attributed, e.g. to the autoregressive residual error component in the model.The simple parameterization of the model error term with three unknown parameters was selected by performing initial fits with different parameterizations, and using sensible initial values refined by a maximum likelihood optimization.The models were then diagnosed by studying the residuals by using normal probability plots and plots of the estimated autocorrelation function.When a good candidate model was found, an MCMC analysis was used to study the uncertainty and identifiability of the variance parameters.Lastly, the in-teresting trend features of the time series were studied by plotting the estimated background level with uncertainty confidence bounds, and drawing realizations from the posterior distribution of the level term and checking the statistical significance of hypothesized features.For example, statistical significance for the trend was assessed by checking whether 95 % posterior probability region included zero value or not.

Estimating trends
A trend is a change in the statistical properties of background state of the system (Chandler and Scott, 2011), the simplest case being linear trend, where, when applicable, we only need to specify the trend coefficient and its uncertainty.companion paper (Kyrölä et al., 2013) the trend and a change in the trend is studied by using a piecewise linear model with a predefined change point.Natural systems evolve continuously in time and it is not always appropriate to approximate the background evolution with a constant, or piecewise linear, trend.Within the state space dynamic linear model framework, the trend can be defined as the change in the estimated background level, i.e. the change in µ t defined in Eq. (3).Posterior sampling from the background level provides an efficient method for studying uncertainties in different trend estimates.
Temporal changes in the system can be studied by visually inspecting the background level and its estimated uncertainty.We can draw samples from the posterior distribution of the level µ t to assess hypotheses about the evolution of the process.For example, in Sect. 3 we study the change in the mean ozone level in 10-year periods.We take into account the uncertainty in the model prediction and in the estimated variance parameters by sampling possible background levels µ t from its posterior uncertainty distribution.We will do this consecutively and, for each sampled realization, calculate the 10-year change in the mean ozone for each time t as trend(t) = µ t+60 − µ t−60 (time units in months).This sample of trends provides us direct way to analyse trends by calculating, for example, the mean 10-year trend with 95 % uncertainty limits.The general procedure is the following: 1. Using MCMC with the Kalman filter likelihood, produce a sample from the marginal posterior probability distribution of the auxiliary parameters defining the error covariance matrix W and model matrix G.
2. Draw one realization of the matrices G and W using the posterior distribution provided by MCMC in the previous step.
3. Simulate one realization of the model state x 1:N using the Kalman simulation smoother assuming fixed G and W from the previous step and calculate trend-related statistics of interest from this realization.
4. Repeat from step 2 to calculate summaries from the posterior distribution of the quantity of interest.

Results
The model parameters were fitted separately to each data set, i.e. to each height interval and zonal band.We performed the analyses using vertical average profile data with both the original interpolated 1 km altitude grid and by forming averaged ozone densities for three altitude regions: 25-35, 35-45, and 45-55 km.The 10 • wide zonal bands start from 60 • S and go to 60 • N. By considering each zonal band independently and summing several altitudes, we have tried to reduce the model to a minimum one that still shows interesting long-term changes and is consistent with our assumptions.Initially, a multivariate estimation was considered by fitting several altitudes and zonal bands together, but this complicated the analyses considerably and did not gain additional insight.In principle, we could use the observations in more refined resolution and model several time series in one estimation step, and even use the individual satellite retrievals instead of spatio-temporal averages.Figure 2 shows an example of our modelling results in the combined altitude-latitude region, 40-50 • N, 35-45 km.The original data are displayed together with DLM estimates and with the estimated mean level component µ t that is used to make statistical inference about the trend.The fits are obtained by a Kalman smoother (Eqs.A1-A7 in Appendix A) and the 95 % uncertainty regions by combination of a simulation smoother (Eqs.A8-A12) and MCMC.A separate panel (on the lower left side) displays the decadal trend obtained from the level component using MCMC analysis to account for the uncertainty in the variance parameters.Looking at the observations in a 10-year perspective, the trend has been statistically significantly negative up to the year 1997, as the grey area stays below zero.After 1997, the 10-year trend does not statistically differ from zero.After the DLM decomposition, the model residual term is assumed to be uncorrelated Gaussian noise.The two lower right panels in Fig. 2 show residual diagnostics.These are used to look for deviations from the modelling assumptions.If the scaled residuals pass the check of being independent and Gaussian with unit variance, then the fit agrees with our assumptions and modelling results are consistent with the data.In this case the residuals are almost perfectly uncorrelated Gaussian with unit variance.The normal probability plot compares empirical quantiles (i.e.order statistics) of the residual to those obtained from the theoretical Gaussian distribution and in the case of normally distributed residuals, the points should stay mostly on a straight line, as they do here.
The results from MCMC analysis for the same data set as in Fig. 2, are shown in Fig. 1.The figure has the prior probability distributions together with marginal posterior distributions as MCMC chain histograms for the estimated auxiliary model parameters (referred as θ in Appendix A).We see that the parameter posterior distributions are mostly determined by the data and not by the prior.The standard deviation of the trend, σ trend , is estimated to be relatively small, which supports the search for smooth background variability.This behaviour is typical for all the fits done.The variability in the seasonal component has more uncertainty, and it also varies more between the fits for different altitude-latitude regions.The autoregressive parameters σ AR and ρ have narrow posterior distributions relative to the prior, i.e. these values can be accurately obtained from the data.
Figures 3 and 4 show the fitted model components for two individual altitudes, 20 and 34 km, both at 10-0 • S. For the seasonal cycle, we plot the sum of the observed seasonal components, u t,1 +u t,2 in Eq. ( 12); for proxies we plot β i z t,i , i.e. the proxy coefficient times the value of the proxy.These fits are provided by the Kalman smoother as in the case of the mean level µ t .At 20 km, we see a strong effect of ENSO.The total variability explained by ENSO is about 20 % (i.e. its effect is between −10 and 10 %), which is similar to the effects of the seasonal variation and of the QBO.Although ENSO was not used by Kyrölä et al. (2013), the fit by linear regression, shown with a solid green line in the upper panel, does agree reasonably well with the DLM level.At 34 km, there is a clear difference between the DLM fit and the piecewise linear analysis with change point fixed at year 1997.Non-linear changes in the background level are not captured by the two-piece linear model.Among the proxy variables, the two QBOs, combined in the figure, have the most significant effect, explaining approximately 10 % of the total variability in ozone at that region.percentages with respect to the total variability in ozone.We see, e.g., that the ENSO index significantly affects the ozone variability only at 20 km near the equator, where the effect is about 20%, also shown in the lowest panel of Fig. 3.The maximal solar effect is always less than 10%.Note that the colour scales differ.

26
Figure 5.The effects of seasonal variation and of the proxy index variables on the variability of ozone.Here 100 % means the observed standard deviation of ozone at each latitude altitude bin.The colouring corresponds to the range of variability during the whole observation period 1984-2011 of each model component given as percentages with respect to the total variability in ozone.We see, e.g. that the ENSO index significantly affects the ozone variability only at 20 km near the equator, where the effect is about 20 %, also shown in the lowest panel of Fig. 3.The maximal solar effect is always less than 10 %.Note that the colour scales differ.
Figure 5 shows the effects of the fitted variable for the seasonal variation and of the proxy index variables on the total observed variability of ozone.We calculated, using the data with 1 km altitude spacing, the range of variability during the whole observation period 1984-2011 of each model component given as percentages with respect to the mean ozone value.The ENSO index has a significant effect (about 20 % of the variability) only at 20 km near the equator.Near 60 km and at higher latitudes, almost all of the ozone variability can be attributed to seasonality.The maximal solar effect is always less than 10 %.
The key results of this paper can be found in Figs. 6 and 7 where we show fits to the ozone time series for 12 latitude belts and for three altitude regions (25-35 km, 35-45 and 45-55 km).We have shown the data points from the merged SAGE II-GOMOS time series, the DLM fit and the slowly varying background level µ t .Overall, it is easy to see that the fits usually follow very accurately the data points.There are a few scattered outliers, but it is difficult to find any specific pattern for these deviations.As in Fig. 2, the solid trend line is the estimated mean level µ t and it is shown with 95 % uncertainty region.The results for the background curve can be roughly divided into three classes.The simplest case, a continuous decay of ozone during the period 1984-2011, is evident in the southernmost latitude band (50-60 • S) in all altitude regions.A continuous decay is also present in 50-60 • N in the lowest altitude region (25-35 km).The second class, a recovery between an initial decay and a final decay, is seen from the equator to 50N in the lowest altitude region (25-35 km) and at 40-60 • N and 35-45 km.The third class covering the rest of the cases shows decay and recovery.In the Southern Hemisphere all altitudes and latitudes (10-50 • S) belong to this class.In the Northern Hemisphere all latitudes in 45-55 km and latitudes from the equator to 40 • N belong to this class.
In Fig. 8 the results of decadal trend analyses are collected and plotted separately for the three altitude ranges.The calculations are explained in Sect.2.5 and as final step 7 in Appendix A. The colouring shows the trend, i.e. the average change in ±5 years around each time axis point.Blue colour means a negative and red a positive trend.The change is most visible at higher altitudes (45-55 and 35-45 km).At the

Figure 9.
To compare DLM results to linear regression analysis in Kyrölä et al. (2013) we reproduce its Figure 16 by using the DLM approach.The left panel shows the change in ozone trends in % per decade between the periods 1997-2011 and 1984-1997.The shaded area indicate regions where the trend difference is nonzero with over 95 % posterior probability.The right panel shows the differences (in % decade −1 ) between the numbers that produce the contour plot on the left panel and the corresponding one in Kyrölä et al. (2013) (Fig. 16).lower altitudes (25-35 km) there are some visible changes, but they are mostly masked by a larger variability and the noise level of the satellite observations at these altitudes.The northernmost (50-60 • N) and the southernmost (60-50 • S zones differ significantly from the other zones.Looking at the individual plots (Figs. 6 and 7), we see some evidence of instrument-related differences and one possible explanation comes from the different non-uniform sampling character-istics, both temporal and spatial, of SAGE II and GOMOS (Toohey et al., 2013).
One particular feature of interest in the data is the suggested stratospheric ozone recovery due to the prohibition of CFC compounds by the Montreal treaty in 1987.Several studies indicate a possible turning point around the year 1997 (Newchurch et al., 2003;Jones et al., 2009;Steinbrecht et al., 2009).According to our analyses, significant changes in the average background level of ozone can be seen at  This suggests that a model with fixed change point in all regions, as used by Kyrölä et al. (2013), does not capture all features in the trends.In some regions there are signs of reduction of ozone after the year 2007.Near the equator the long-term background changes are mostly masked by irregular variations in the ozone concentration.Kyrölä et al. (2013) analysed the trend change point by calculating several fixed piecewise linear regression analyses and choosing the one that provided largest difference between before and after the point change in the linear trend.In the DLM approach the trend can change its value continuously and we can analyse directly where the most likely change points are.One use of a DLM model would be the identification of possible trend change points, and provide hypotheses on which a traditional static regression analysis would be performed.The results of the piecewise linear trend analysis in Kyrölä et al. (2013) are in agreement to some extent with the results here, but the DLM analysis suggests that the trend change point depends on latitude and altitude, as shown in Fig. 8. Also, there is some evidence against using a simple two-piece linear model to describe the background changes.

Comparison with piecewise linear analysis
To further compare the methods, we studied the difference in trend between the years from 1984 to 1997 and from 1997 to 2011 using the DLM method.Figure 16 of Kyrölä et al. (2013) corresponds to our Fig. 9 (left panel).Overall, the DLM method produces somewhat larger differences (i.e.faster ozone recovery) in before-and after-1997 trends.This difference varies and is at most 3 % / 10 years.An opposite effect is seen in the southern most latitude band (50-60 • S), where the DLM method suggests stronger decline of ozone than the linear method.The uncertainty (as a standard deviation of the estimated change) in DLM trend difference is about 1.5 % / 10 year on the average.For the piecewise linear model it is approximately 0.5 % larger.Because of this relatively large uncertainty, the differences between DLM and piecewise linear trend analysis are mostly not statistically significant.The right panel in Fig. 9 shows the difference in the results between the linear and the DLM approach.The analysis for the before-and-after-1997 trend differs mostly at the lowest and highest latitudes and at the equator.At the southernmost zone the difference is the largest and can be attributed to the ability of the DLM method to dynamically adapt to changing seasonal patterns in the observations, see Fig. 7.However, this might be caused by the different spatiotemporal sampling of the two instruments.Further analysis is beyond the scope of this paper -see the discussion in the companion paper (Kyrölä et al., 2013).
Figure 10 shows six examples of the DLM approach vs. the piecewise linear model fits from Kyrölä et al. (2013).In five of the cases the conclusions disagree in some respect.In these cases the DLM trend has more variability than the fixed-point model trend -the linear model may miss some important features in the data, and the rigidness of the piecewise linear model may cause spurious results.In some cases the change point is later than at the beginning of the year 1997.However, at those regions where the one-change-point model is valid, the conclusions are the same with both modelling approaches.

Conclusions
We have shown that a dynamic linear model (DLM) is well suited for modelling ozone time series.In contrast to some classical time series analyses, DLM does not require stationarity, it allows for missing observations and takes uncertainties in the observations into account.By using Markov chain Monte Carlo (MCMC) simulation analysis, the uncertainty in the structural variance parameters can be accounted for.The state space method directly includes a model error term, which makes the analysis more robust to mis-specification of the model.The analysis allows full statistical uncertainty quantification, and it is extendible to more refined analyses, if those seem necessary.Here, a relatively straightforward and conceptually simple analysis reveals interesting features and, also, validates the more ad hoc choices used in piecewise linear global analyses, such as in the companion paper by Kyrölä et al. (2013).
Trend analysis can be a delicate matter and it is always challenging to give causal explanations.With a properly set up DLM model we can detect smooth changes in the background state.By using proxy variables we can filter out the effect of known external forcing, such as the solar effect.The DLM analysis provides a method to detect and quantify trends, but the statistical model itself does not provide explanations.It can verify that the observations are consistent with the selected model.Model diagnostics will eventually falsify wrong models and other badly selected prior specifications.
We used a local linear trend model with two harmonic functions for the seasonal effect and four proxy time series for the solar flux, Quasi-Biennial Oscillations and ENSO.An autoregressive model component was used to account for possible residual autocorrelation.The results show a statistically significant change point in the combined SAGE II-GOMOS time series approximately after year 1997 at altitudes 35-55 km and mid-latitudes, between 50-20 • S and 20-50 • N.This change point in time varies with latitude and altitude.There are locations where the estimated changes are opposite to those expected and the length of the data set is still short relative to some cycles of natural variability in the atmospheric processes.At the southernmost zones analysed here (60-50 • S) the DLM method shows larger decline of ozone.This might be attributed to the sampling differences of the two instruments, which is better handled by the DLM method.In many regions the behaviour of the mean ozone concentration is too complicated to allow for a simple piecewise linear description.Compared to the two-piece fixed change point linear model used in the companion paper, we see stronger recovery of ozone in many regions -the difference is up to 3 % / decade.parameter, with the state x 1:N uncertainty accounted for, i.e. integrated out.Fortunately, this likelihood is available as a by-product of the Kalman filter calculations for each fixed θ.By the assumed Markov property of the state space equations, this marginal likelihood function can be evaluated sequentially as a product of the individual time wise marginal likelihoods as p(y 1:N |θ) = p(y 1 |θ) N t=2 p(y t |y 1:t−1 , θ).
(A13) where x t is the one-step-ahead mean prediction of the state and C y,t is the covariance matrix of predicted observation, both obtained by the Kalman filter formulas, given above in Eqs.(A1) and (A3).
Step 5 With this likelihood function and with a prior distribution for the parameter vector θ, a MCMC simulation can be performed to produce a sample of parameter values from the marginal posterior distribution p(θ|y 1:N ).Thus, we can estimate and study the uncertainty in θ by the MCMC.We use an adaptive version of the random walk Metropolis-Hastings (MH) algorithm by Haario et al. (2006).It is a sequential Monte Carlo algorithm that proposes a value for the parameter vector from a suitably constructed proposal distribution and then either accepts or rejects it depending on the ratio of the posterior distribution at the proposed value to that calculated using the previously accepted value.By the Bayes formula, the posterior distribution is, up to a normalizing constant, equal to the likelihood times the prior.The basic idea behind MCMC is the fact that the acceptance probability depends only on the ratio of the posterior at two parameter values, and this will cancel out the hard-to-calculate proportional constant in the Bayes formula for the posterior distribution.See Gamerman (2006) for more details on the MCMC methodology.A MCMC Matlab toolbox used in the calculations is available from http://helios.fmi.fi/~lainema/mcmc/.
Step 6 As final steps in our DLM calculations, we can combine the simulation smoother and the MCMC outputs to produce samples from the joint posterior distribution p(x 1:N , θ|y 1:N ), as well as from the marginal posterior distribution of the state, p(x 1:N |y 1:N ), where the uncertainty in θ has been accounted for.
Step 7 Lastly, a sample from this last distribution is used for the trend analyses by calculating trend statistics that depend on the estimated model states, as explained in Sect.2.5.

Fig. 1 .Figure 1 .
Fig. 1. Results of MCMC analysis for the four unknown model parameter: trend change standard deviation σ trend , seasonal change standard deviation σ trend , autoregressive standard deviation σ AR , and the autoregressive coefficient ρ.The histograms are estimates of the corresponding marginal posterior probability densities obtained by MCMC analysis with 10 000 simulation steps.The parameters are for 40 • N-50 • N, 35-45 km asin Fig.2.The prior probability densities given by dashed lines are log-normal for the standard deviations and truncated Gaussian for the AR parameter ρ.For the analysis the ozone observation have been scaled by their observed standard deviation.This means that the values given in the x axes are also scaled so that the overall variance of the ozone observations is one.

Figure 2 .
Figure 2. DLM fit for average ozone at 40-50 • N, 35-45 km.In the upper panel the dots are the observations used in the analysis, the solid line following the observations is the DLM fit obtained by a Kalman filter.The smooth solid line is the background level component of the model with 95 % probability envelope.The beginning of GOMOS observations as well as the end of SAGE II are shown by red vertical lines.In the lower left panel an analysis of the 10-year trend is performed, showing the mean estimated trend and its 95 % probability envelope.The lower right panel shows model diagnostic analyses on the residuals by an estimated autocorrelation function (ACF) and by a normal probability plot.In the autocorrelation function plot, dashed horizontal lines correspond to the approximative region where the coefficients do not significantly (95 %) deviate from zero.The unit of lag is months.

Fig. 3 .Figure 3 .
Fig. 3. Some of the DLM model components for 10 • S-0 • S, 20 km.The first panel has the observations and the level components as in Fig. 2. In addition, the piecewise linear regression model used by Kyrölä et al. (2013) is drawn with a solid green line.The second panel downwards shows model residuals, scaled by the estimated residual standard deviation, i.e. the variability in the y axis direction should be approximately standard Gaussian N (0, 1).The other panels show means and 95% probability envelopes (the grey shading) of the DLM components related to seasonality and to the proxy variables, solar, Quasi-Biennial Oscillation, and ENSO, i.e. the original proxy variables multiplied by the estimated regression coefficients.For the proxy variables, the y axis scaling is the percentage of the total ozone variability, i.e. the ozone observations were scaled by their standard deviation, individually for each data set analysed.At 20 km we see a large number of missing observations, some possible outliers, and a strong effect of the ENSO index.There is only a slight difference between DLM and linear models, however, although the latter does not include ENSO.

Fig. 4 .Figure 4 .
Fig. 4. Some of the DLM model components for 10 • S-0 • S, 34 km, see Fig.3for explanation.We see noticeable difference between DLM and linear regression fits drawn with a solid green line on the top panel.

Fig. 5 .
Fig. 5.The effects of seasonal variation and of the proxy index variables on the variability of ozone.Here 100% means the observed standard deviation of ozone at each latitude altitude bin.The colouring corresponds to the range of variability during the whole observation period 1984-2011 of each model component given as

Fig. 6 .Figure 6 .Fig. 7 .
Fig.6.All fits collected for the northern hemisphere.The three altitude regions as columns, northern hemisphere 10 • zonal bands as rows starting from the northernmost zone.The smooth solid curve shows the estimated background level with 95% probability envelope.The dots are the observations used in the analysis.The solid line following the observations is the DLM fit obtained by Kalman filter formulas.

Figure 7 .
Figure 7.All fits collected for the Southern Hemisphere.The three altitude regions as columns, Southern Hemisphere 10 • zonal bands as rows starting from the equatorial zone.See Fig. 6 for explanation.

Fig. 8 .Figure 8 .Fig. 9 .
Fig. 8. Ten-year trend in percentage change / year from the DLM fits for each altitude region and zonal band.One individual trend analysis is shown in the lower left panel of Fig. 2.

Fig. 10 .Figure 10 .
Fig. 10.Comparing DLM and the linear model.Six examples of DLM vs. piecewise linear model fits using data with 1 km altitude resolution.The solid black line is the piecewise linear trend from Kyrölä et al. (2013).The smooth solid blue line is the background level component of the DLM model with 95% probability envelope.At 40 • N-50 • N and 44 km the methods agree.In the other panels we see differences of various degrees.

Figure 11 .
Figure 11.A flowchart showing the dependencies of the DLM calculations explained in Appendix A.
The rationale behind G seas(k) is that if we know the harmonic u t,k = a k cos(k2π/12 t) + b k sin(k2π/12 t) and, as an auxiliary state, its conjugate u * t,k = −a k sin(k2π/12 t)+b k cos(k2π/12 t), with some constants a k and b k , we can update the state with

Table 1 .
Specification of priors for auxiliary model parameters estimated by MCMC.The prior distributions for the model error standard deviation parameters were log-normal, logN (log(µ), σ 2 ) with values of µ and σ given in the table.For simplicity, the same relative priors are used in all the altitude-latitude regions.In Fig.1the prior distributions are plotted together with MCMC chain histograms that estimate the posterior distributions for one altitude-latitude region.
These individual one-step-ahead predictive distributions of the observations are provided by the Kalman filter.In the case of a linear Gaussian model this likelihood is proportional to (ignoring a constant that does not depend on the model parameters) − F t x t ) T C −1 y,t (y t − F t x t ) + log( C y,t ) ,