ARMA Forecasting for Non-Gaussian Time-Series Data Using Copulas | by Sarem Seitz | Jun, 2022
It’s all fun and games until your time-series is not Normal anymore
ARMA (AutoRegressive — Moving Average) models are arguably the most popular approach to time-series forecasting. Unfortunately, plain ARMA is made for Gaussian distributed data only. On the one hand, you can often still use ARMA by transforming the raw data. On the other hand, this typically makes probabilistic forecasts quite tedious.
One approach to apply ARMA to non-Normal data are Copula models. Roughly, the latter allow us to exchange the Gaussian marginal for any other continuous distribution. At the same time, they preserve the implicit time-dependency between observations that is imposed by ARMA.
If this sounds confusing, I suggest reading the next paragraph carefully. Also, you might want to read some external sources for a deeper understanding, too.
Informally, Copulas (or Copulae if you are a Latin hardliner) define joint cumulative distribution functions (c.d.f.) for unit-uniform random variables. Formally, we can describe this as
That property alone is quite unspectacular as uniform random variables are not very expressive for practical problems. However, an important result in probability theory will make things more interesting.
The probability integral transform states that we can transform any continuous random variable to a uniform one by plugging it into its own c.d.f.:
We can verify this empirically for a standard Normal example:
As the inverse of a c.d.f. is the quantile function, we can easily invert this transformation. Even cooler, we can transform a uniform random variable to any continuous random variable via
This inverse transformation will become relevant later on.
Combining Copulas and the inverse probability transform
In conjunction with Copulas, this allows us to separate the marginal distributions from the dependency structure of joint random variables.
A concrete example: Consider two random variables, X and Y with standard Gamma and Beta marginal distributions, i.e.
With the help of a Copula and the probability integral transform, we can now define a joint c.d.f over both variables such that we preserve their marginal distributions:
Introducing the Gaussian Copula
So far, we haven’t specified any Copula function yet. A simplistic one is the Gaussian Copula, which is defined as follows:
If we combine this with the Gamma-Beta example from before, we get the following Gaussian Copula joint c.d.f.:
The implicit rationale behind this approach can be described in three steps:
- Transform the Gamma and Beta marginals into Uniform marginals via the respective c.d.f.s
- Transform the Uniform marginals into standard Normal marginals via the quantile functions
- Define the joint distribution via the multivariate Gaussian c.d.f. with zero mean, unit variance and non-zero covariance (covariance matrix
R)
By inverting these steps, we can easily sample from a bi-variate random variable that has the above properties. I.e. standard Gamma/Beta marginals with Gaussian Copula dependencies:
- Draw a sample from a bi-variate Gaussian with mean zero, unit variance and non-zero covariance (covariance matrix
R). You now have two correlated standard Gaussian variables. - Transform both variables with the standard Gaussian c.d.f. — you now have two correlated Uniform variables. (= probability integral transform)
- Transform these variables with the standard Beta and Gamma quantile functions — you now have a pair of correlated GammaBeta variables. (= inverse probability integral transform)
Notice that we could drop the zero-mean, unit-variance assumption on the multivariate Gaussian. In that case we would have to adjust the Gaussian c.d.f. to the corresponding marginals in order to keep the integral probability transform valid.
Since we are only interested in the dependency structure (i.e. covariances), standard Gaussian marginals are sufficient and easier to deal with.
Now let us sample some data in Julia:
Congratulations, you have just sampled from your first Copula model!
But wait — I want to fit a model!
Let’s say we observed the above data without the underlying generating process. We only presume that we know that Gamma-Beta marginals and a Gaussian copula are a good choice. How could we fit the model parameters (i.e. ‘learn’ them in the Machine Learning world)?
As often for statistical models, Maximum Likelihood is a good approach. However, we need a density function for that, so what do we do? We already found out, that a Copula model describes a valid c.d.f. for continuous marginals. Thus, we can derive the corresponding probability density by taking derivatives:
If you carefully look at the last line, you see that the joint density equals the product of the marginal densities times the Copula density. Recall that the product of the marginal densities alone would imply independency. Thus, the Copula density serves as a correction factor to account for the actual dependency between the marginals.
This structure allows us to separate the marginal distributions from the their dependency structure.
Now, for the Gaussian Copula, one can prove the following Copula density function:
In order to estimate the model parameters for the Gamma-Beta example, we simply need to plug this function into the full density. Then, we take logarithms and get the target log-likelihood function for optimization.
Finally, we can return to our initial problem. For this example, we will focus on the stationary ARMA(1,1) model:
For a time-series with T observations, we can derive the unconditional, stationary distribution (see e.g. here):
Informally, the unconditional distribution considers a fixed-length time-series as a single, multivariate random vector. As a consequence, it doesn’t matter whether we are sampling from the unconditional distribution or the usual ARMA equations (for an equally long time-series) themselves.
In some instances, such as this one, the unconditional distribution is easier to work with.
Also, notice that the unconditional marginal distributions (the distributions of the y_t's) are the same regardless of the time-lag we are looking at. In fact, we have zero-mean Gaussians with variance equal to the auto-covariance function at zero.
Next, let us define:
The transformed covariance matrix now implies unit variance while preserving the dependency structure of the unconditional time-series. Literally, we have just derived the correlation matrix but let us stick to the idea of a standardized covariance matrix.
If we plug this back into a Gaussian copula — let us call the ARMA(1,1) Copula. Now, we could use the ARMA(1,1) Copula dependency structure together with any continuous marginal distribution. For example, we could define
i.e. the unconditional marginals are Exponential-distributed with rate parameter 0.5. Putting everything together, we obtain the following unconditional density:
Let us combine everything so far and plot an example:
Clearly, the samples from the Copula model are not Gaussian anymore. In fact, we observe a single draw from a ARMA(1,1) Copula with Exponential-distributed marginals.
Parameter estimation with Maximum Likelihood
So far, we have only been able to simulate a time-series from the ARMA(1,1) Copula model. In order to fit the model, we will apply Maximum Likelihood. When using Copulas for cross-sectional data, it is usually possible to separate fitting the marginal distributions from fitting the Copula. Unfortunately, this does not work here.
As we only observe one realization of the process per marginal, fitting a distribution based on the marginals alone is impossible. Rather, we now need to optimize both the marginals and the copula at once. This begs the additional difficulty of having to deal with the marginal’s parameters inside the marginal’s c.d.f..
Namely, our Maximum likelihood objective looks as follows:
Optimizing this can become quite ugly as derivatives with respect to a c.d.f.’s parameters are usually fairly complex. Luckily, the Exponential distribution is quite simple and respective derivatives are easily found. Even better, the Optim.jl package can optimize our log-likelihood via finite differences without requiring any derivatives at all.
If we chose another distribution than the Exponential, finite differences might not suffice. In that case, we would have to either implement the c.d.f. derivatives by hand or hope that ChainRules.jl can handle them for us.
Also, we transform our model parameters to the correct domains via exp and tanh instead of applying Box constraints in the Optim optimizer. This worked reasonably accurate and fast here:
Now, let us evaluate the result. For the Exponential distribution, the estimated parameter should be close to the true parameter. Regarding the latent ARMA parameters, we primarily need the estimated auto-covariance to be close to ground-truth. This is indeed the case here:
Forecasting with the Copula model
Finally, we want to use our model to produce actual forecasts. Due to the Copula construction, we can derive the conditional forecast density in closed form. As we will see however, mean and quantile forecasts need to be calculated numerically.
First, recall how the Copula model defines a joint density over all ‘training’-observations:
In order to forecast a conditional density at h steps ahead, we simply need to follow standard probability laws:
This boils down to the ratio of two Copula evaluations times the marginal density evaluated at the target point. However, we still need to find a way to use this equation to calculate a mean forecast and a forecast interval.
As the density is arguably fairly complex, we won’t even try to derive any of these values in closed form. Rather, we use numerical methods to find the target quantities.
For the mean, we simply use quadrature to approximate the usual integral:
with U a sufficiently large value to capture most of the probability mass (approximation up to infinity is obviously not possible).
For the forecast interval, we use the 90% prediction interval. Thus, we need to find the 5% and the 95% quantiles of the conditional density. This can be done via another approximation, this time through an Ordinary Differential Equation:
For a derivation of this formula, see, for example, here. Integrating the ODE from zero up to the target quantile yields the respective target quantile value. The latter can be done numerically via DifferentialEquations.jl.
With this, we can finally calculate the forecast and plot the result:
The looks indeed quite reasonable and the forecast appears to converge to a stable distribution as we predict further ahead into the future.
As we have seen, Copulas make it possible to extend well-known models to non-Gaussian data. This allowed us to transfer the simplicity of the ARMA model to Exponential marginals that were only defined for positive values.
One complication arises when the size of the observed time-series becomes very large. In that case, the unconditional covariance matrix will scale poorly and the model fitting step will likely become impossible.
Then, we need to find a computationally more efficient solution. One possible approach are Implicit Copulas which define a Copula density through a chain of conditional densities.
Of course, there are many other ways to integrate Copulas into classical statistical and Machine Learning models. For the latter, research is still a little sparse. However, I strongly believe that there is at least some potential for a modern application of these classic statistical objects.
[1] Hamilton, James Douglas. Time series analysis. Princeton university press, 2020.
[2] Nelsen, Roger B. An introduction to copulas. Springer Science & Business Media, 2007.
[3] Smith, Michael Stanley. Implicit copulas: An overview. Econometrics and Statistics, 2021.
It’s all fun and games until your time-series is not Normal anymore
ARMA (AutoRegressive — Moving Average) models are arguably the most popular approach to time-series forecasting. Unfortunately, plain ARMA is made for Gaussian distributed data only. On the one hand, you can often still use ARMA by transforming the raw data. On the other hand, this typically makes probabilistic forecasts quite tedious.
One approach to apply ARMA to non-Normal data are Copula models. Roughly, the latter allow us to exchange the Gaussian marginal for any other continuous distribution. At the same time, they preserve the implicit time-dependency between observations that is imposed by ARMA.
If this sounds confusing, I suggest reading the next paragraph carefully. Also, you might want to read some external sources for a deeper understanding, too.
Informally, Copulas (or Copulae if you are a Latin hardliner) define joint cumulative distribution functions (c.d.f.) for unit-uniform random variables. Formally, we can describe this as
That property alone is quite unspectacular as uniform random variables are not very expressive for practical problems. However, an important result in probability theory will make things more interesting.
The probability integral transform states that we can transform any continuous random variable to a uniform one by plugging it into its own c.d.f.:
We can verify this empirically for a standard Normal example:
As the inverse of a c.d.f. is the quantile function, we can easily invert this transformation. Even cooler, we can transform a uniform random variable to any continuous random variable via
This inverse transformation will become relevant later on.
Combining Copulas and the inverse probability transform
In conjunction with Copulas, this allows us to separate the marginal distributions from the dependency structure of joint random variables.
A concrete example: Consider two random variables, X and Y with standard Gamma and Beta marginal distributions, i.e.
With the help of a Copula and the probability integral transform, we can now define a joint c.d.f over both variables such that we preserve their marginal distributions:
Introducing the Gaussian Copula
So far, we haven’t specified any Copula function yet. A simplistic one is the Gaussian Copula, which is defined as follows:
If we combine this with the Gamma-Beta example from before, we get the following Gaussian Copula joint c.d.f.:
The implicit rationale behind this approach can be described in three steps:
- Transform the Gamma and Beta marginals into Uniform marginals via the respective c.d.f.s
- Transform the Uniform marginals into standard Normal marginals via the quantile functions
- Define the joint distribution via the multivariate Gaussian c.d.f. with zero mean, unit variance and non-zero covariance (covariance matrix
R)
By inverting these steps, we can easily sample from a bi-variate random variable that has the above properties. I.e. standard Gamma/Beta marginals with Gaussian Copula dependencies:
- Draw a sample from a bi-variate Gaussian with mean zero, unit variance and non-zero covariance (covariance matrix
R). You now have two correlated standard Gaussian variables. - Transform both variables with the standard Gaussian c.d.f. — you now have two correlated Uniform variables. (= probability integral transform)
- Transform these variables with the standard Beta and Gamma quantile functions — you now have a pair of correlated GammaBeta variables. (= inverse probability integral transform)
Notice that we could drop the zero-mean, unit-variance assumption on the multivariate Gaussian. In that case we would have to adjust the Gaussian c.d.f. to the corresponding marginals in order to keep the integral probability transform valid.
Since we are only interested in the dependency structure (i.e. covariances), standard Gaussian marginals are sufficient and easier to deal with.
Now let us sample some data in Julia:
Congratulations, you have just sampled from your first Copula model!
But wait — I want to fit a model!
Let’s say we observed the above data without the underlying generating process. We only presume that we know that Gamma-Beta marginals and a Gaussian copula are a good choice. How could we fit the model parameters (i.e. ‘learn’ them in the Machine Learning world)?
As often for statistical models, Maximum Likelihood is a good approach. However, we need a density function for that, so what do we do? We already found out, that a Copula model describes a valid c.d.f. for continuous marginals. Thus, we can derive the corresponding probability density by taking derivatives:
If you carefully look at the last line, you see that the joint density equals the product of the marginal densities times the Copula density. Recall that the product of the marginal densities alone would imply independency. Thus, the Copula density serves as a correction factor to account for the actual dependency between the marginals.
This structure allows us to separate the marginal distributions from the their dependency structure.
Now, for the Gaussian Copula, one can prove the following Copula density function:
In order to estimate the model parameters for the Gamma-Beta example, we simply need to plug this function into the full density. Then, we take logarithms and get the target log-likelihood function for optimization.
Finally, we can return to our initial problem. For this example, we will focus on the stationary ARMA(1,1) model:
For a time-series with T observations, we can derive the unconditional, stationary distribution (see e.g. here):
Informally, the unconditional distribution considers a fixed-length time-series as a single, multivariate random vector. As a consequence, it doesn’t matter whether we are sampling from the unconditional distribution or the usual ARMA equations (for an equally long time-series) themselves.
In some instances, such as this one, the unconditional distribution is easier to work with.
Also, notice that the unconditional marginal distributions (the distributions of the y_t's) are the same regardless of the time-lag we are looking at. In fact, we have zero-mean Gaussians with variance equal to the auto-covariance function at zero.
Next, let us define:
The transformed covariance matrix now implies unit variance while preserving the dependency structure of the unconditional time-series. Literally, we have just derived the correlation matrix but let us stick to the idea of a standardized covariance matrix.
If we plug this back into a Gaussian copula — let us call the ARMA(1,1) Copula. Now, we could use the ARMA(1,1) Copula dependency structure together with any continuous marginal distribution. For example, we could define
i.e. the unconditional marginals are Exponential-distributed with rate parameter 0.5. Putting everything together, we obtain the following unconditional density:
Let us combine everything so far and plot an example:
Clearly, the samples from the Copula model are not Gaussian anymore. In fact, we observe a single draw from a ARMA(1,1) Copula with Exponential-distributed marginals.
Parameter estimation with Maximum Likelihood
So far, we have only been able to simulate a time-series from the ARMA(1,1) Copula model. In order to fit the model, we will apply Maximum Likelihood. When using Copulas for cross-sectional data, it is usually possible to separate fitting the marginal distributions from fitting the Copula. Unfortunately, this does not work here.
As we only observe one realization of the process per marginal, fitting a distribution based on the marginals alone is impossible. Rather, we now need to optimize both the marginals and the copula at once. This begs the additional difficulty of having to deal with the marginal’s parameters inside the marginal’s c.d.f..
Namely, our Maximum likelihood objective looks as follows:
Optimizing this can become quite ugly as derivatives with respect to a c.d.f.’s parameters are usually fairly complex. Luckily, the Exponential distribution is quite simple and respective derivatives are easily found. Even better, the Optim.jl package can optimize our log-likelihood via finite differences without requiring any derivatives at all.
If we chose another distribution than the Exponential, finite differences might not suffice. In that case, we would have to either implement the c.d.f. derivatives by hand or hope that ChainRules.jl can handle them for us.
Also, we transform our model parameters to the correct domains via exp and tanh instead of applying Box constraints in the Optim optimizer. This worked reasonably accurate and fast here:
Now, let us evaluate the result. For the Exponential distribution, the estimated parameter should be close to the true parameter. Regarding the latent ARMA parameters, we primarily need the estimated auto-covariance to be close to ground-truth. This is indeed the case here:
Forecasting with the Copula model
Finally, we want to use our model to produce actual forecasts. Due to the Copula construction, we can derive the conditional forecast density in closed form. As we will see however, mean and quantile forecasts need to be calculated numerically.
First, recall how the Copula model defines a joint density over all ‘training’-observations:
In order to forecast a conditional density at h steps ahead, we simply need to follow standard probability laws:
This boils down to the ratio of two Copula evaluations times the marginal density evaluated at the target point. However, we still need to find a way to use this equation to calculate a mean forecast and a forecast interval.
As the density is arguably fairly complex, we won’t even try to derive any of these values in closed form. Rather, we use numerical methods to find the target quantities.
For the mean, we simply use quadrature to approximate the usual integral:
with U a sufficiently large value to capture most of the probability mass (approximation up to infinity is obviously not possible).
For the forecast interval, we use the 90% prediction interval. Thus, we need to find the 5% and the 95% quantiles of the conditional density. This can be done via another approximation, this time through an Ordinary Differential Equation:
For a derivation of this formula, see, for example, here. Integrating the ODE from zero up to the target quantile yields the respective target quantile value. The latter can be done numerically via DifferentialEquations.jl.
With this, we can finally calculate the forecast and plot the result:
The looks indeed quite reasonable and the forecast appears to converge to a stable distribution as we predict further ahead into the future.
As we have seen, Copulas make it possible to extend well-known models to non-Gaussian data. This allowed us to transfer the simplicity of the ARMA model to Exponential marginals that were only defined for positive values.
One complication arises when the size of the observed time-series becomes very large. In that case, the unconditional covariance matrix will scale poorly and the model fitting step will likely become impossible.
Then, we need to find a computationally more efficient solution. One possible approach are Implicit Copulas which define a Copula density through a chain of conditional densities.
Of course, there are many other ways to integrate Copulas into classical statistical and Machine Learning models. For the latter, research is still a little sparse. However, I strongly believe that there is at least some potential for a modern application of these classic statistical objects.
[1] Hamilton, James Douglas. Time series analysis. Princeton university press, 2020.
[2] Nelsen, Roger B. An introduction to copulas. Springer Science & Business Media, 2007.
[3] Smith, Michael Stanley. Implicit copulas: An overview. Econometrics and Statistics, 2021.
Denial of responsibility! Techno Blender is an automatic aggregator of the all world’s media. In each content, the hyperlink to the primary source is specified. All trademarks belong to their rightful owners, all materials to their authors. If you are the owner of the content and do not want us to publish your materials, please contact us by email – [email protected]. The content will be deleted within 24 hours.