If the variance of the observations in the series is not constant across the range of observations it may be useful to apply a variance-stabilizing transformation to the series. A common situation is for the variance to increase with the magnitude of the observations and in this case typical transformations used are the log or square root transformation. A range-mean plot or standard deviation-mean plot provides a quick and easy way of detecting non-constant variance and of choosing, if required, a suitable transformation. These are plots of either the range or standard deviation of successive groups of observations against their means.

These may be used to simplify the structure of a time series.

First-order differencing, i.e., forming the new series will remove a linear trend. First-order seasonal differencing eliminates a fixed seasonal pattern.

These operations reflect the fact that it is often appropriate to model a time series in terms of changes from one value to another. Differencing is also therefore appropriate when the series has something of the nature of a random walk, which is by definition the accumulation of independent changes.

Differencing may be applied repeatedly to a series, giving where d and D are the orders of differencing. The derived series wt will be shorter, of length N=n-d-s×D, and extend for t=1+d+s×D,…,n.

Given that a series has (possibly as a result of simplifying by differencing operations) a homogeneous appearance throughout its length, fluctuating with approximately constant variance about an overall mean level, it is appropriate to assume that its statistical properties are stationary. For most purposes the correlations ρk between terms xt,xt+k or wt,wt+k separated by lag k give an adequate description of the statistical structure and are estimated by the sample autocorrelation function (ACF) rk, for k=1,2,….

The information in the autocorrelations, ρk, may be presented in a different light by deriving from them the coefficients of the partial autocorrelation function (PACF) ϕk,k, for k=1,2,…. ϕk,k which measures the correlation between xt and xt+k conditional upon the intermediate values xt+1,xt+2,…,xt+k-1. The corresponding sample values ϕ^k,k give further assistance in the selection of ARIMA models.

Both autocorrelation function (ACF) and PACF may be rapidly computed, particularly in comparison with the time taken to estimate ARIMA models.

The partial autocorrelation coefficient ϕk,k is determined as the final parameter in the minimum variance predictor of xt in terms of xt-1,xt-2,…,xt-k, where ek,t is the prediction error, and the first subscript k of ϕk,i and ek,t emphasises the fact that the parameters will alter as k increases. Moderately good estimates ϕ^k,i of ϕk,i are obtained from the sample autocorrelation function (ACF), and after calculating the partial autocorrelation function (PACF) up to lag L, the successive values v1,v2,…,vL of the prediction error variance estimates, vk=varek,t, are available, together with the final values of the coefficients ϕ^k,1,ϕ^k,2,…,ϕ^k,L. If xt has nonzero mean, x-, it is adequate to use xt-x- in place of xt in the prediction equation.

In theory, the parameters of an ARIMA model are determined by a sufficient number of autocorrelations ρ1,ρ2,…. Using the sample values r1,r2,… in their place it is usually (but not always) possible to solve for the corresponding ARIMA parameters.

These are rapidly computed but are not fully efficient estimates, particularly if moving-average parameters are present. They do provide useful preliminary values for an efficient but relatively slow iterative method of estimation. This is based on the least-squares principle by which parameters are chosen to minimize the sum of squares of the innovations at, which are regenerated from the data using (2), or the reverse of (3) and (4) in the case of seasonal models.

Lack of knowledge of terms on the right-hand side of (2), when t=1,2,…,maxp,q, is overcome by introducing q unknown series values w0,w1,…,wq-1 which are estimated as nuisance parameters, and using correction for transient errors due to the autoregressive terms. If the data w1,w2,…,wN=w is viewed as a single sample from a multivariate Normal density whose covariance matrix V is a function of the ARIMA model parameters, then the exact likelihood of the parameters is The least-squares criterion as outlined above is equivalent to using the quadratic form as an objective function to be minimized. Neglecting the term -12logV yields estimates which differ very little from the exact likelihood except in small samples, or in seasonal models with a small number of whole seasons contained in the data. In these cases bias in moving-average parameters may cause them to stick at the boundary of their constraint region, resulting in failure of the estimation method.

Approximate standard errors of the parameter estimates and the correlations between them are available after estimation.

The model residuals, a^t, are the innovations resulting from the estimation and are usually examined for the presence of autocorrelation as a check on the adequacy of the model.

An ARIMA model is particularly suited to extrapolation of a time series. The model equations are simply used for t=n+1,n+2,… replacing the unknown future values of at by zero. This produces future values of wt, and if differencing has been used this process is reversed (the so-called integration part of ARIMA models) to construct future values of xt.

Forecast error limits are easily deduced.

This process requires knowledge only of the model orders and parameters together with a limited set of the terms at-i,et-i,wt-i,xt-i which appear on the right-hand side of the models (3) and (4) (and the differencing equations) when t=n. It does not require knowledge of the whole series.

We call this the state set. It is conveniently constituted after model estimation. Moreover, if new observations xn+1,xn+2,… come to hand, then the model equations can easily be used to update the state set before constructing forecasts from the end of the new observations. This is particularly useful when forecasts are constructed on a regular basis. The new innovations an+1,an+2,… may be compared with the residual standard deviation, σa, of the model used for forecasting, as a check that the model is continuing to forecast adequately.

Exponential smoothing is a relatively simple method of short term forecasting for a time series. A variety of different smoothing methods are possible, including; single exponential, Brown's double exponential, linear Holt (also called double exponential smoothing in some references), additive Holt–Winters and multiplicative Holt–Winters. The choice of smoothing method used depends on the characteristics of the time series. If the mean of the series is only slowly changing then single exponential smoothing may be suitable. If there is a trend in the time series, which itself may be slowly changing, then linear Holt smoothing may be suitable. If there is a seasonal component to the time series, e.g., daily or monthly data, then one of the two Holt–Winters methods may be suitable.

For a time series yt, for t=1,2,…,n, the five smoothing functions are defined by the following:

• Single Exponential Smoothing

  mt = α yt + 1-α mt-1 y^t+f = mt var y^t+f = varεt 1+ f-1 α2

• Brown Double Exponential Smoothing

  mt = α yt + 1-α mt-1 rt = α mt - mt-1 + 1-α rt-1 y^ t+f = mt + f-1 + 1 / α rt var y^ t+f = varεt 1+ ∑ i=0 f-1 2α+ i-1 α2 2

• Linear Holt Smoothing

  mt = α yt + 1-α mt-1 + ϕ rt-1 rt = γ mt - mt-1 + 1-γ ϕ rt-1 y^ t+f = mt + ∑ i=1 f ϕi rt var y^ t+f = var εt 1+ ∑ i=1 f-1 α + α γ ϕ ϕi-1 ϕ-1 2

• Additive Holt–Winters Smoothing

  mt = α yt - s t-p + 1-α m t-1 +ϕ r t-1 rt = γ mt - m t-1 + 1-γ ϕ rt-1 st = β yt - mt + 1-β s t-p y^ t+f = mt + ∑ i=1 f ϕi rt + s t-p var y^ t+f = var εt 1+ ∑ i=1 f-1 ψi2 ψi = 0 if ​i≥f α + αγϕ ϕ i -1 ϕ-1 if ​ i mod p≠0 α + α γ ϕ ϕi -1 ϕ-1 +β 1-α otherwise

• Multiplicative Holt–Winters Smoothing

  mt = α yt / s t-p + 1-α m t-1 +ϕ r t-1 rt = γ mt - m t-1 + 1-γ ϕ r t-1 st = β yt / mt + 1-β s t-p y^ t+f = mt + ∑ i=1 f ϕi rt × s t-p var y^ t+f = var εt ∑ i=0 ∞ ∑ j=0 p-1 ψ j+ip s t+f s t+f-j 2

  and ψ is defined as in the additive Holt–Winters smoothing,

where mt is the mean, rt is the trend and st is the seasonal component at time t with p being the seasonal order. The f-step ahead forecasts are given by y^t+f and their variances by var y^ t+f . The term var εt  is estimated as the mean deviation.

The parameters, α, β and γ control the amount of smoothing. The nearer these parameters are to one, the greater the emphasis on the current data point. Generally these parameters take values in the range 0.1 to 0.3. The linear Holt and two Holt–Winters smoothers include an additional parameter, ϕ, which acts as a trend dampener. For 0.0<ϕ<1.0 the trend is dampened and for ϕ>1.0 the forecast function has an exponential trend, ϕ=0.0 removes the trend term from the forecast function and ϕ=1.0 does not dampen the trend.

For all methods, values for α, β, γ and ψ can be chosen by trying different values and then visually comparing the results by plotting the fitted values along side the original data. Alternatively, for single exponential smoothing a suitable value for α can be obtained by fitting an ARIMA0,1,1 model. For Brown's double exponential smoothing and linear Holt smoothing with no dampening, (i.e., ϕ=1.0), suitable values for α and, in the case of linear Holt smoothing, γ can be obtained by fitting an ARIMA0,2,2 model. Similarly, the linear Holt method, with ϕ≠1.0, can be expressed as an ARIMA1,2,2 model and the additive Holt–Winters, with no dampening, (ϕ=1.0), can be expressed as a seasonal ARIMA model with order p of the form ARIMA0,1,p+10,1,0. There is no similar procedure for obtaining parameter values for the multiplicative Holt–Winters method, or the additive Holt–Winters method with ϕ≠1.0. In these cases parameters could be selected by minimizing a measure of fit using nonlinear optimization.

For a series x1,x2,…,xn this is defined as the scaling factor now being chosen in order that i.e., the spectrum indicates how the sample variance (σx2) of the series is distributed over components in the frequency range 0≤ω≤π.

It may be demonstrated that f*ω is equivalently defined in terms of the sample ACF rk of the series as where ck=σx2rk are the sample autocovariance coefficients.

If the series xt does contain a deterministic sinusoidal component of amplitude R, this will be revealed in the sample spectrum as a sharp peak of approximate width π/n and height n/2πR2. This is called the discrete part of the spectrum, the variance R2 associated with this component being in effect concentrated at a single frequency.

If the series xt has no deterministic components, i.e., is purely stochastic being stationary with autocorrelation function (ACF) rk, then with increasing sample size the expected value of f*ω converges to the theoretical spectrum – the continuous part where γk are the theoretical autocovariances.

The sample spectrum does not however converge to this value but at each frequency point fluctuates about the theoretical spectrum with an exponential distribution, being independent at frequencies separated by an interval of 2π/n or more. Various devices are therefore employed to smooth the sample spectrum and reduce its variability. Much of the strength of spectral analysis derives from the fact that the error limits are multiplicative so that features may still show up as significant in a part of the spectrum which has a generally low level, whereas they are completely masked by other components in the original series. The spectrum can help to distinguish deterministic cyclical components from the stochastic quasi-cycle components which produce a broader peak in the spectrum. (The deterministic components can be removed by regression and the remaining part represented by an ARIMA model.)

A large discrete component in a spectrum can distort the continuous part over a large frequency range surrounding the corresponding peak. This may be alleviated at the cost of slightly broadening the peak by tapering a portion of the data at each end of the series with weights which decay smoothly to zero. It is usual to correct for the mean of the series and for any linear trend by simple regression, since they would similarly distort the spectrum.

The estimate is calculated directly from the sample autocovariances ck as the smoothing being induced by the lag window weights wk which extend up to a truncation lag M which is generally much less than n. The smaller the value of M, the greater the degree of smoothing, the spectrum estimates being independent only at a wider frequency separation indicated by the bandwidth b which is proportional to 1/M. It is wise, however, to calculate the spectrum at intervals appreciably less than this. Although greater smoothing narrows the error limits, it can also distort the spectrum, particularly by flattening peaks and filling in troughs.

The unsmoothed sample spectrum is calculated for a fine division of frequencies, then averaged over intervals centred on each frequency point for which the smoothed spectrum is required. This is usually at a coarser frequency division. The bandwidth corresponds to the width of the averaging interval.

An important tool for investigating how an input series xt affects an output series yt is the sample cross-correlation function (CCF) rxyk, for k=0,1,2,…, between the series. If xt and yt are (jointly) stationary time series this is an estimator of the theoretical quantity The sequence ryxk, for k=0,1,2,…, is distinct from rxyk, though it is possible to interpret When the series yt and xt are believed to be related by a transfer function (TF) model, the CCF is determined by the impulse response function (IRF) v0,v1,v2,… and the autocorrelation function (ACF) of the input xt.

In the particular case when xt is an uncorrelated series or white noise (and is uncorrelated with any other inputs): and the sample CCF can provide an estimate of vk: where sy and sx are the sample standard deviations of yt and xt, respectively.

In theory the IRF coefficients vb,…,vb+p+q determine the parameters in the TF model, and using v~k to estimate v~k it is possible to solve for preliminary estimates of δ1,δ2,…,δp, ω0,ω1,…,ωq.

In general an input series xt is not white noise, but may be represented by an ARIMA model with innovations or residuals at which are white noise. If precisely the same operations by which at is generated from xt are applied to the output yt to produce a series bt, then the transfer function relationship between yt and xt is preserved between bt and at. It is then possible to estimate The procedure of generating at from xt (and bt from yt) is called prewhitening or filtering by an ARIMA model. Although at is necessarily white noise, this is not generally true of bt.

A multi-input model may be used to forecast the output series provided future values (possibly forecasts) of the input series are supplied.

Construction of the forecasts requires knowledge only of the model orders and parameters, together with a limited set of the most recent variables which appear in the model equations. This is called the state set. It is conveniently constituted after model estimation. Moreover, if new observations yn+1,yn+2,… of the output series and xn+1,xn+2,… of (all) the independent input series become available, then the model equations can easily be used to update the state set before constructing forecasts from the end of the new observations. The new innovations an+1,an+2,… generated in this updating may be used to monitor the continuing adequacy of the model.

In many time series applications it is desired to calculate the response (or output) of a transfer function (TF) model for a given input series.

Smoothing, detrending, and seasonal adjustment are typical applications. You must specify the orders and parameters of a TF model for the purpose being considered. This may then be applied to the input series.

Again, problems may arise due to ignorance of the input series values prior to the observation period. The transient errors which can arise from this may be substantially reduced by using ‘backforecasts’ of these unknown observations.

As in the case of a univariate time series, it may be useful to simplify the series by differencing operations which may be used to remove linear or seasonal trends, thus ensuring that the resulting series to be used in the model estimation is stationary. It may also be necessary to apply transformations to the individual components of the multivariate series in order to stabilize the variance. Commonly used transformations are the log and square root transformations.

Multivariate analogues of the autocorrelation and partial autocorrelation functions are available for analysing a set of k time series, xi,1,xi,2,…,xi,n, for i=1,2,…,k, thereby making it possible to obtain some understanding of a suitable VARMA model for the observed series.

It is assumed that the time series have been differenced if necessary, and that they are jointly stationary. The lagged correlations between all possible pairs of series, i.e., are then taken to provide an adequate description of the statistical relationships between the series. These quantities are estimated by sample auto- and cross-correlations rijl. For each l these may be viewed as elements of a (lagged) autocorrelation matrix.

Thus consider the vector process xt (with elements xit) and lagged autocovariance matrices Γl with elements of σiσjρijl where σi2=varxi,t. Correspondingly, Γl is estimated by the matrix Cl with elements sisjrijl where si2 is the sample variance of xit.

For a series with short-term cross-correlation only, i.e., rijl is not significant beyond some low lag q, then the pure vector MAq model, with no autoregressive parameters, i.e., p=0, is appropriate.

The cross-correlation structure of a stationary multivariate time series may often be represented by a model with a small number of parameters belonging to the VARMA class. If the stationary series wt has been derived by transforming and/or differencing the original series xt, then wt is said to follow the VARMA model: where εt is a vector of uncorrelated residual series (white noise) with zero mean and constant covariance matrix Σ, ϕ1,ϕ2,…,ϕp are the p autoregressive (AR) parameter matrices and θ1,θ2,…,θq are the q moving-average (MA) parameter matrices. If wt has a nonzero mean μ, then this can be allowed for by replacing wt,wt-1,… by wt-μ,wt-1-μ,… in the model.

A series generated by this model will only be stationary provided restrictions are placed on ϕ1,ϕ2,…,ϕp to avoid unstable growth of wt. These are stationarity constraints. The series εt may also be usefully interpreted as the linear innovations in wt, i.e., the error if wt were to be predicted using the information in all past values wt-1,wt-2,…, provided also that θ1,θ2,…,θq satisfy what are known as invertibility constraints. This allows the series εt to be generated by rewriting the model equation as The method of maximum likelihood (ML) may be used to estimate the parameters of a specified VARMA model from the observed multivariate time series together with their standard errors and correlations.

The residuals from the model may be examined for the presence of autocorrelations as a check on the adequacy of the fitted model.

Forecasts of the series may be constructed using a multivariate version of the univariate method. Efficient methods are available for updating the forecasts each time new observations become available.

This is defined by It may be demonstrated that this is equivalently defined in terms of the sample cross-correlation function (CCF), rxyk, of the series as where cxyk=sxsyrxyk is the cross-covariance function.

The cross-spectrum is specified by its real part or cospectrum cf*ω and imaginary part or quadrature spectrum qf*ω, but for the purpose of interpretation the cross-amplitude spectrum and phase spectrum are useful:

If the series xt and yt contain deterministic sinusoidal components of amplitudes Ry,Rx and phases ϕy,ϕx at frequency ω, then A*ω will have a peak of approximate width π/n and height n/2πRyRx at that frequency, with corresponding phase ϕ*ω=ϕy-ϕx. This supplies no information that cannot be obtained from the two series separately. The statistical relationship between the series is better revealed when the series are purely stochastic and jointly stationary, in which case the expected value of fxy*ω converges with increasing sample size to the theoretical cross-spectrum where γxyk=covxt,yt+k. The sample spectrum, as in the univariate case, does not converge to the theoretical spectrum without some form of smoothing which either implicitly (using a lag window) or explicitly (using a frequency window) averages the sample spectrum fxyω* over wider bands of frequency to obtain a smoothed estimate f^xyω.

If there is no statistical relationship between the series at a given frequency, then fxyω=0, and the smoothed estimate f^xyω, will be close to 0. This is assessed by the squared coherency between the series: where f^xxω is the corresponding smoothed univariate spectrum estimate for xt, and similarly for yt. The coherency can be treated as a squared multiple correlation. It is similarly invariant in theory not only to simple scaling of xt and yt, but also to filtering of the two series, and provides a useful test statistic for the relationship between autocorrelated series. Note that without smoothing, so the coherency is 1 at all frequencies, just as a correlation is 1 for a sample of size 1. Thus smoothing is essential for cross-spectrum analysis.

The estimate of the cross-spectrum is calculated from the sample cross-variances as The lag window wk extends up to a truncation lag M as in the univariate case, but its centre is shifted by an alignment lag S usually chosen to coincide with the peak cross-correlation. This is equivalent to an alignment of the series for peak cross-correlation at lag 0, and reduces bias in the phase estimation.

The selection of the truncation lag M, which fixes the bandwidth of the estimate, is based on the same criteria as for univariate series, and the same choice of M and window shape should be used as in univariate spectrum estimation to obtain valid estimates of the coherency, gain, etc., and test statistics.

The computations are exactly as for smoothing of the univariate spectrum except that allowance is made for an implicit alignment shift S between the series.

To improve the stability of the computations the square root algorithm is used. One recursion of the square root covariance filter algorithm which can be summarized as follows: where U is an orthogonal transformation triangularizing the left-hand pre-array to produce the right-hand post-array, St is the lower triangular Cholesky factor of the state covariance matrix Pt+1∣t, Qt1/2 and Rt1/2 are the lower triangular Cholesky factor of the covariance matrices Q and R and H1/2 is the lower triangular Cholesky factor of the covariance matrix of the residuals. The relationship between the Kalman gain matrix, Kt, and Gt is given by To improve the efficiency of the computations when the matrices At,Bt and Ct do not vary with time the system can be transformed to give a simpler structure. The transformed state vector is U*X where U* is the transformation that reduces the matrix pair A,C to lower observer Hessenberg form. That is, the matrix U* is computed such that the compound matrix is a lower trapezoidal matrix. The transformations need only be computed once at the start of a series, and the covariance matrices Qt and Rt can still be time-varying.

If the state space model contains unknown parameters, θ, these can be estimated using maximum likelihood (ML). Assuming that Wt and Vt are normal variates the log-likelihood for observations Yt, for t=1,2,…,n, is given by Optimal estimates for the unknown model parameters θ can then be obtained by using a suitable optimizer method to maximize the likelihood function.

Once the model has been fitted forecasting can be performed by using the one-step-ahead prediction equations. The one-step-ahead prediction equations can also be used to ‘jump over’ any missing values in the series.

Many commonly used time series models can be written as state space models. A univariate ARMAp,q model can be cast into the following state space form: where r=maxp,q+1.

The representation for a k-variate ARMAp,q series (VARMA) is very similar to that given above, except now the state vector is of length kr and the ϕ and θ are now k×k matrices and the 1s in A, B and C are now the identity matrix of order k. If p<r or q+1<r then the appropriate ϕ or θ matrices are set to zero, respectively.

Since the compound matrix is already in lower observer Hessenberg form (i.e., it is lower trapezoidal with zeros in the top right-hand triangle) the invariant Kalman filter algorithm can be used directly without the need to generate a transformation matrix U*.

Rather than modelling the mean (for example using regression models) or the autocorrelation (by using ARMA models) there are circumstances in which the variance of a time series needs to be modelled. This is common in financial data modelling where the variance (or standard deviation) is known as volatility. The ability to forecast volatility is a vital part in deciding the risk attached to financial decisions like portfolio selection. The basic model for relating the variance at time t to the variance at previous times is the autoregressive conditional heteroskedastic (ARCH) model. The standard ARCH model is defined as where ψt is the information up to time t and ht is the conditional variance.