Linear Stochastic Models in Discrete and Continuous Time

Abstract

The econometric data to which autoregressive moving-average models are commonly applied are liable to contain elements from a limited range of frequencies. If the data do not cover the full Nyquist frequency range of $[0,\pi ]$ radians, then severe biases can occur in estimating their parameters. The recourse should be to reconstitute the underlying continuous data trajectory and to resample it at an appropriate lesser rate. The trajectory can be derived by associating sinc fuction kernels to the data points. This suggests a model for the underlying processes. The paper describes frequency-limited linear stochastic differential equations that conform to such a model, and it compares them with equations of a model that is assumed to be driven by a white-noise process of unbounded frequencies. The means of estimating models of both varieties are described.

1. Introduction

The paper considers the means by which a linear stochastic differential equation (LSDE) can be inferred from a series of observations taken at finite intervals in time. The objective is to explore the means of estimating models that are driven by a frequency-limited white noise and to compare them with the conventional models that are assumed to be driven by a continuous white noise comprising the increments of a Wiener process. Such a process has a spectrum that is unbounded in frequency.

The paper represents a sequel to a previous paper (Pollock 2012), which described the hazards of estimating a discrete-time autoregressive moving-average (ARMA) model on the basis of data that have been sampled at a rate that exceeds the minimum rate that is needed in order to capture the information contained in the underlying continuous process. It was shown that the resulting estimates can suffer from significant biases.

To obtain consistent estimates of the ARMA parameters, it is necessary to re-sample the data at a rate that is attuned to the maximum frequency of the data, so as to equate that frequency to the Nyquist frequency of $\pi$ radians per sampling interval. The re-sampling can be achieved by reconstituting a continuous process from the sampled ordinates.

According to the Shannon–Nyquist sampling theorem, this can be done by attaching a scaled sinc function kernel to each of the sampled ordinates. The theorem provides a model for a frequency-limited stochastic processes in continuous time, which has been explored in author’s previous paper. It is now intended to represent such a model by a stochastic differential equation.

The Shannon–Nyquist theorem presupposes a doubly-infinite domain, which is necessary in order to support the sinc functions, which extend indefinitely in both directions. In the case of a finite sample, the sinc function must be replaced by a circularly wrapped sinc function, or Dirichlet kernel, which is applied to the data by a circular convolution.

In the present paper, it is shown that a Dirichlet kernel interpolation is equivalent to a Fourier synthesis, involving a finite sum of continuous trigonometrical functions. In the case of a doubly-infinite sequence, the sum transmutes into an integral, and Fourier synthesis ceases to be practicable. Then, one must have recourse to the Shannon–Nyquist formulation.

A white-noise process that consists of the increments of a Wiener process can be depicted as the limiting case of a frequency-limited process that arises when the sampling interval tends to zero. As the interval diminishes, the sinc functions become increasingly concentrated on the sample points. In the limit, they become impulses that can be represented by Dirac’s delta function. In the process, the sample ordinates that scale the sinc functions become infinitesimal.

It will be observed that a sinc function is differentiable to any order. The frequency-limited process, which consists of a superposition of sinc functions, is likewise differentiable. In the limiting case, where the sampling interval tends to zero, the trajectory retains its continuity, but it ceases to be differentiable. Despite this difficulty, many authors continue to employ the derivatives operator D in their formulations of stochastic differential equations.

As Bergstrom (1984) has observed, in the course of a rigorous analysis of these matters, this is allowable provided that it is understood that the process thereby denoted is defined only by the properties of its integral. We shall exploit this allowance, if only for notational convenience. It allows us to treat the frequency-limited white noise and the white noise of unbounded frequencies together in the same notation. Much of the embarrassment can be overcome by adopting the spectral representation of the white-noise process that drives the system.

Bergstrom has been the progenitor of a distinct branch of continuous-time analysis that has been pursued by econometricians, of which a brief account should be given.

A major preoccupation of the econometricians is with structural multivariate continuous-time models, in which the equations are liable to represent the aggregate behaviour of economic agents at any point in time.

To this end, it is important to observe the distinction between stock variables, such money stock, exchange rates, and asset prices, which are measured at points in time, and flow variables, such as consumers’ expenditure, income and aggregate investment, which are measured as integrals of a rate of flow, taken over the observation period.

Bergstrom, whose primary concern was with vector first-order continuous-time processes, showed that the inclusion of flow variables in the model induces a moving-average term.

Bergstrom’s method of estimation was to derive a so-called exact discrete linear model (EDLM), which corresponds in discrete time to the continuous-time system, and from which the parameters of the latter could be inferred. He was able to overcome the problem of the non-observability of the derivatives in a second-order system by a process of substitution, which eliminates them from the system.

Chambers (1999) has extended the methodology of Bergstrom, to encompass differential equations of higher orders. He has overcome the problem of the non-observability of the derivatives of the variables by extending Bergstrom’s method of substitution.

In the process, he has shown that a continuous-time model in stock variables with an autoregressive order of p and a moving-average order of $q<p$, which is an LSDE$(p,q)$, will have a discrete-time counterpart of orders p and $p-1$, which is an EDLM$(p,p-1)$. He has also shown that, if the continuous system includes flow variables, then the counterpart discrete-time model will be an EDLM($p,p$). For a detailed derivation, see also McCrorie (2000).

In this paper, attention is confined to univariate LSDE$(p,q)$ models in stock variables. The principal concern of the paper is with the case of a continuous-time model in which the maximum frequency is the Nyquist frequency of $\pi$ radian per sample interval.

This model, which has a one-to-one correspondence with a discrete-time ARMA model, will be described as a continuous-time ARMA (CARMA) model. The acronym CARMA will be reserved exclusively for such models, despite the fact that it is commonly used to denote an LSDE model driven by a white-noise process of unbounded frequencies.

The correspondence of the LDSE$(p,q)$ model and EDLM$(p,p-1)$ model, which is an ARMA$(p,p-1)$ model of a special provenance, is demonstrated, in this paper, in the context of the partial-fraction decomposition of the LSDE.

This renders it as a set first-order components that are operating in parallel. The components have real-valued parameters, or they have complex-valued parameters, which come in conjugate pairs. The method of eliminating the non-observable derivatives by a process of substitution is thereby avoided, albeit that it continues to be the requisite method in the case of a multivariate process.

It is also evident, from the partial-fraction decomposition of an improper rational function, which includes a quotient term, that an LSDE$(p,p)$ model, which might accommodate flow variables, will have a discrete-time EDLM$(p,p)$ counterpart.

The plan of the paper is as follows. In Section 2, the means of reconstructing a continuous signal from its sampled ordinates is considered, both in the case of a finite sequence and in the case of a doubly-infinite sequence, which is the subject of the sampling theorem of Nyquist (1924, 1928) and of Shannon (1949). Section 3 provides the essential results concerning discrete-time ARMA models and continuous-time LSDE models in which the forcing function is a white-noise process comprising the increments of a Wiener process, and Section 4 concerns the autocovariance functions and the spectral density functions of these models.

A central result of the paper is contained in Section 5, which derives a version of the frequency-limited ARMA process in continuous-time. This is achieved by attaching sinc functions to the ordinates of the discrete-time ARMA process. However, it is also required to derive a frequency-limited LSDE to represent the process that generates the continuous trajectory of the signal.

The principle of impulse invariance serves to establish an immediate one-to-one relationship between the parameters of the ARMA model and those of the frequency-limited LSDE model. To find the parameters of an LSDE that is unbounded in frequency, the alternative principle of autocovariance equivalence must be invoked. This involves an iterative procedure, which is described in Section 6.

This section, which is followed by a brief summary of the paper, also provides examples of LSDE models in which the forcing function consists of the increments of a Wiener process, and of continuous-time CARMA models that are driven by frequency-limited white noise.

Section 2 and Section 3 rely on the Fourier transforms of discrete and continuous functions. An account of these, which includes the matters of sampling and aliasing, which can affect functions in both the time domain and the frequency domain, has been provided by Pollock (2018). in an appendix to the paper. A more extensive treatment of Fourier analysis is given by Pollock (1999).

2. Reconstrusting and Resampling the Signal

Consider a continuous real-valued signal $x\left(t\right)$ defined on the interval $[0,T)$, and imagine that it is sampled at unit intervals at the points $t=0,1,\dots ,T-1$. Then, the discrete Fourier transform of the sequence $x\left(t\right)=\{x_t;t=0,1,\dots ,T-1\}$ delivers a sequence of complex-valued coefficients $\xi \left(j\right)=\{{\xi }_j;j=0,1,\dots T-1\}$; and there is a one-to-one correspondence between the two sequences, which is expressed as

$$\left\{x_t=\sum _{j=0}^{T-1}{\xi }_j{e}^{\mathrm{i}{\omega }_{j}t}\right\}⟷\left\{{\xi }_j=\frac{1}{T}\sum _{t=0}^{T-1}{x}_t{e}^{-\mathrm{i}{\omega }_{j}t}\right\},\mathrm{with}{\omega }_j=\frac{2\pi j}{T}.$$

(1)

The expression on the LHS is described as the spectral representation of the sequence, or as its Fourier synthesis, and the expression on the RHS defines its Fourier analysis.

Here, ${\omega }_j$ takes T equally spaced values in the interval $[0,2\pi )$. However, since $exp\left(\mathrm{i}\omega \right)=cos\left(\omega \right)+\mathrm{i}sin\left(\omega \right)$ is a $2\pi$-periodic function, it is possible to replace the positive frequencies $(2\pi -\omega )\in [\pi ,2\pi )$ by the negative frequencies $-\omega \in [-\pi ,0)$. Then, since $cos\left(\omega \right)=\{exp(\mathrm{i}\omega )+exp(-\mathrm{i}\omega \left)\right\}/2$ and $sin\left(\omega \right)=-\mathrm{i}\{exp(\mathrm{i}\omega )-exp(-\mathrm{i}\omega \left)\right\}/2$, the data sequence can be expressed as

$$x\left(t\right)=\sum _{j=0}^{T-1}{\xi }_j{e}^{\mathrm{i}{\omega }_{j}t}=\sum _{j=0}^{[T/2]}\{{\alpha }_{j}cos\left({\omega }_{j}t\right)+{\beta }_{j}sin\left({\omega }_{j}t\right)\},$$

(2)

where ${\alpha }_j={\xi }_j+{\xi }_{T-j}$ and ${\beta }_j=\mathrm{i}({\xi }_j-{\xi }_{T-j})$ are real-valued coefficients and where $exp\left\{{\omega }_{T-j}\right\}=exp\{-{\omega }_j\}$. Furthermore, $[T/2]$ stands for the integral part of $T/2$, which enables the expression to accommodate an odd-valued sample size.

By replacing in (2) the index $t=0,1,\dots ,T-1$ of the sample by a continuous variable $t\in [0,T)$, a trajectory is generated that interpolates the sampled elements according to the method of Fourier interpolation.

When T is even-valued, the nonzero frequencies of the trigonometric functions of (2) range from ${\omega }_1=2\pi /T$ to ${\omega }_{T/2}=\pi$. The former corresponds to one cycle in the time spanned by the data, and the latter is the Nyquist frequency, which is the frequency of the function ${(-1)}^t$ and which represents the highest frequency that is observable within sampled data.

If, in its trigonometrical representation, the underlying signal contains elements at frequencies in excess of $\pi$, then these frequencies will be confounded in the sampled sequence with frequencies falling within the observable Nyquist range of $[0,\pi ]$ radians per sampling interval, which are described as their aliases.

If the limiting frequency within the signal is $\pi$ radians per unit interval, then it should be possible, via Equation (2), to reconstruct the trajectory of the signal over the interval $[0,T)$ by letting t vary continuously in the interval.

The Fourier synthesis of (2) depicts the segment of the signal over the interval $[0,T)$ as a single cycle of a continuous periodic function. For this to be an accurate representation of the signal, it is necessary that $x_{T-1}\simeq x_0$, which is to say that there must be no significant discontinuity where one replication of the segment ends and another begins.

In effect, the segment representing the transition between $x_{T-1}$ and $x_0$ is required to be synthesised from the available Fourier frequencies. To replicate a jump discontinuity or “saltus” in the underlying signal would require an infinite Fourier sum comprising an unbounded set of frequencies.

Various recourses are available to ensure that this requirement is met. Thus, it may be appropriate to replace the sampled ordinates by their deviations from a line that passes through the points $x_0$ and $x_{T-1}$ or that passes close to them. Then, once it has been generated by the Fourier interpolation of the deviations, the continuous trajectory may be added back to the line. This technique may be employed whenever there is a significant trend in the data.

An alternative recourse, which is available when there is no significant trend and which can also be applied to the detrended data, is to interpolate a segment of pseudo data between $x_{T-1}$ and $x_0$ so as to ensure a continuous transition between these two points. Once a continuous trajectory has been created, the segment based on the pseudo data can be discarded. Both of these recourses are available in the computer program IDEOLOG.PAS, available on the author’s website, and examples have been provided in Pollock (2015).

Using the expression for ${\xi }_j$ from (1), the representation of the continuous signal $x\left(t\right)$ over the interval $[0,T)$ that is based on T elements of sampled data can be written as

$$x\left(t\right)=\sum _{j=0}^{T-1}\left\{\frac{1}{T}\sum _{k=0}^{T-1}{x}_k{e}^{-\mathrm{i}{\omega }_{j}k}\right\}{e}^{\mathrm{i}{\omega }_{j}t}=\frac{1}{T}\sum _{k=0}^{T-1}{x}_k\left\{\sum _{j=0}^{T-1}{e}^{\mathrm{i}{\omega }_j(t-k)}\right\}.$$

(3)

Here, it is understood than ${\omega }_{(T/2)+j}={\omega }_{j-(T/2)}$. Therefore, in the case where T is odd-valued, the inner summation of the final expression gives rise to the following Dirichlet Kernel:

$${\phi }_n^{\circ }\left(t\right)=\frac{1}{T}\sum _{j=-n}^n{e}^{\mathrm{i}{\omega }_{j}t}=\frac{sin\left(\{T/2\}{\omega }_{1}t\right)}{Tsin({\omega }_{1}t/2)}\mathrm{with}n=\frac{T-1}{2}.$$

(4)

In the case where T is even-valued, there is a marginally more complicated expression for the Dirichlet kernel that has been given, for example, by Pollock (2009), where the derivations are provided in both odd and even cases. Thus, the continuous signal of finite duration can also be expressed in terms of the Dirichlet kernel:

$$x\left(t\right)=\sum _{t=0}^{T-1}{x}_k{\phi }^{\circ }(t-k).$$

(5)

This method of recovering the continuous signal from its sampled ordinates is described as kernel interpolation.

Since theoretical time-series analysis commonly presupposes that the data constitute a doubly-infinite sequence, it is also appropriate to consider the limiting form of the discrete Fourier transform as the size of the sample grows indefinitely. This is the discrete-time Fourier transform. Thus, in the limit, Equation (1) becomes

$$\left\{x_t=\frac{1}{2\pi }{\int }_{-\pi }^{\pi }\xi \left(\omega \right)e^{\mathrm{i}\omega t}d\omega \right\}⟷\xi \left(\omega \right)=\sum _{t=-\infty }^{\infty }{x}_t{e}^{-\mathrm{i}\omega t}.$$

(6)

By allowing t to vary continuously, the expression for an underlying continuous signal of an indefinite duration that is limited in frequency to $\pi$ radians becomes

$$x\left(t\right)=\frac{1}{2\pi }{\int }_{-\pi }^{\pi }\left\{\sum _{k=-\infty }^{\infty }{x}_k{e}^{-\mathrm{i}k\omega }\right\}{e}^{\mathrm{i}\omega t}d\omega =\frac{1}{2\pi }\sum _{k=-\infty }^{\infty }{x}_k\left\{{\int }_{-\pi }^{\pi }{e}^{\mathrm{i}\omega (t-k)}\right\}d\omega .$$

(7)

The integral on the RHS is evaluated as

$${\int }_{-\pi }^{\pi }{e}^{\mathrm{i}\omega (t-k)}d\omega =2\frac{sin\left\{\pi \right(t-k\left)\right\}}{t-k}.$$

(8)

Putting this into the RHS gives a weighted sum of sinc functions:

$$x\left(t\right)=\sum _{k=-\infty }^{\infty }{x}_{k}sin\left\{\pi \right(t-k\left)\right\}\pi (t-k)=\sum _{k=-\infty }^{\infty }{x}_k\phi (t-k),$$

(9)

where

$$\phi (t-k)=\frac{sin\left\{\pi \right(t-k\left)\right\}}{\pi (t-k)}$$

(10)

is the sinc function centred on $t=k$. Equation (9) represents sinc function interpolation, and it corresponds to the classic result of Shannon (1949).

It is easy to see that the Dirichlet kernel represents a circularly wrapped version of the sinc function. Compare

$$\phi \left(\omega \right)={\int }_{-\infty }^{\infty }\phi \left(t\right)e^{-\mathrm{i}\omega t}dt⟷\phi \left(t\right)={\int }_{-\pi }^{\pi }{e}^{-\mathrm{i}\omega t}d\omega ,$$

(11)

which relates the sinc function to its Fourier transform, with

$${\phi }^{\circ }\left({\omega }_j\right)={\int }_{-T/2}^{T/2}{\phi }^{\circ }\left(t\right)e^{-\mathrm{i}{\omega }_{j}t}dt⟷{\phi }^{\circ }\left(t\right)=\sum _{t=-n}^n{e}^{-\mathrm{i}{\omega }_{j}t};n=\frac{T-1}{2},$$

(12)

which is the analogous relationship for the Dirichlet kernel. The equality ${\phi }^{\circ }\left({\omega }_j\right)=\phi \left({\omega }_j\right)$ must prevail at each of the Fourier frequencies ${\omega }_j$, and, for the corresponding integrals with respect to t to be equal, it is necessary that

$${\phi }^{\circ }\left(t\right)=\sum _{j=-\infty }^{\infty }\phi (t+jT);t=0,\pm 1,\dots ,\pm T-12.$$

(13)

Thus, the Dirichlet kernel would be created by wrapping the sinc function around a circle of circumference T and by adding the overlying ordinates.

An ARMA process is supported in the frequency domain on the Nyquist interval, which is $(0,\pi ]$, if one is considering trigonometric functions, or $(-\pi ,\pi ]$, if one is considering complex exponential functions. It has been show by Pollock (2012) that severe biases can arise in the estimation of such a model whenever the process is bounded by a value ${\omega }_c<\pi$ less than the Nyquist frequency.

An appropriate recourse, in that case, may be to reconstitute the continuous signal from its sampled ordinates by Fourier interpolation and to resample this at the rate of one observation in each interval of $\pi /{\omega }_c$ time units. In the case where $k=\pi /{\omega }_c$ is an integer, it is appropriate to subsample the data by taking one in every k sample points. The effect, in either case, will be to expand the spectrum of the data to cover the full Nyquist interval.

A different circumstance arises when the maximum frequency ${\omega }_c>\pi$ of the underlying continuous signal exceeds the Nyquist value. The continuous signal and its Fourier transform are related as follows:

$$x\left(t\right)=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }{e}^{\mathrm{i}\omega t}\xi \left(\omega \right)d\omega ⟷\xi \left(\omega \right)={\int }_{-\infty }^{\infty }{e}^{-\mathrm{i}\omega t}x\left(t\right)dt.$$

(14)

For an element of the sampled sequence of $\{x_t;t=0,\pm 1,\pm 2,\dots \}$, there is

$$x_t=\frac{1}{2\pi }{\int }_{-\pi }^{\pi }{e}^{\mathrm{i}\omega t}{\xi }_S\left(\omega \right)d\omega ⟷{\xi }_S\left(\omega \right)=\sum _{k=-\infty }^{\infty }{x}_k{e}^{-ik\omega }.$$

(15)

Therefore, at $x_t=x\left(t\right)$, there is

$$x_t=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }{e}^{\mathrm{i}\omega t}\xi \left(\omega \right)d\omega =\frac{1}{2\pi }{\int }_{-\pi }^{\pi }{e}^{\mathrm{i}\omega t}{\xi }_S\left(\omega \right)d\omega .$$

(16)

The equality of the two integrals implies that

$${\xi }_S\left(\omega \right)=\sum _{j=-\infty }^{\infty }\xi (\omega +2j\pi ).$$

(17)

If the domain of $\xi \left(\omega \right)$ is not limited in frequency to $[-\pi ,\pi )$, then ${\xi }_S\left(\omega \right)\ne \xi \left(\omega \right)$ and aliasing will occur. It can be seen that, in this case, the function ${\xi }_S\left(\omega \right)$ arises from wrapping $\xi \left(\omega \right)$ around a circle in the frequency domain with a circumference of $2\pi$.

This is an essential result in communications engineering, and it can be seen as an adjunct of the Nyquist–Shannon sampling theorem. Econometricians were alerted to it by Granger and Hatanaka (1964) and reminded of it by Granger and Newbold (1977). It was presented to them in the form above by Fishman (1969).

3. Models in Discrete and Continuous Time

A linear stochastic differential equation of orders p and $q<p$, denoted by LSDE$(p,q)$, that is driven by a process $\zeta \left(t\right)$ that is mean-square differentiable, is specified by the equation

$$\begin{array}{ccc}\hfill ({\varphi }_0{D}^p+{\varphi }_1{D}^{p-1}+\cdots +{\varphi }_p)y\left(t\right)& =\hfill & ({\theta }_0{D}^q+{\theta }_1{D}^{q-1}+\cdots +{\theta }_q)\zeta \left(t\right)\hfill \\ \hfill \mathrm{or}\varphi \left(D\right)y\left(t\right)& =\hfill & \theta \left(D\right)\zeta \left(t\right),\hfill \end{array}$$

(18)

wherein D is the derivative operator such that $Dx\left(t\right)=dx\left(t\right)/dt$. Replacing D within the rational operator $\theta \left(D\right)/\varphi \left(D\right)$ by the complex variable s gives rise to the so-called transfer function $\theta \left(s\right)/\varphi \left(s\right)$, which is the Laplace transform of the operator. See, for example, James (2004, p. 168).

An appropriate specification for $\zeta \left(t\right)$ is a frequency-limited white-noise process of the sort to be considered later in Section 5, albeit that the limiting frequency may be virtually unbounded. However, it is also necessary to consider the case where $\zeta \left(t\right)$ is a continuous white-noise process comprising the increments of a Wiener process that are unbounded in frequency and which is therefore non-differentiable.

Such a process has spectral representation that can be denoted by

$$\zeta \left(t\right)={\int }_{-\infty }^{\infty }{e}^{\mathrm{i}\omega t}dZ\left(\omega \right),$$

(19)

where $Z\left(\omega \right)$ is a process of a zero mean and of uncorrelated complex-valued increments with the property that

$$E\left\{dZ\left(t\right)dZ(t-\tau )\right\}={\sigma }_{\zeta }^2\delta \left(\tau \right),$$

(20)

where $\delta \left(\tau \right)$ is Dirac’s delta. The value of $dZ\left(\omega \right)$ can be regarded as a constant in any realisation. Therefore, the derivatives operator D may be applied, under the integral, to the complex exponential term. This justifies us, for example, in writing

$$y\left(t\right)=\theta \left(D\right)\zeta \left(t\right)={\int }_{-\infty }^{\infty }\theta \left(\mathrm{i}\omega \right)e^{\mathrm{i}\omega t}dZ\left(\omega \right),$$

(21)

where $\theta \left(\mathrm{i}\omega \right)={\theta }_0{\left(\mathrm{i}\omega \right)}^q+{\theta }_1{\left(\mathrm{i}\omega \right)}^{q-1}+\cdots +{\theta }_q$. The result extends to the case of the rational operator $\theta \left(D\right)/\varphi \left(D\right)$ by virtue of its series expansion. The analysis of the LSDE process can proceed on the understudying that, whenever $\zeta \left(t\right)$ is to be regarded as a white-noise process of unbounded frequencies, such substitutions are implicit.

It will be observed that the function

$$\psi \left(\mathrm{i}\omega \right)=\frac{\theta \left(\mathrm{i}\omega \right)}{\varphi \left(\mathrm{i}\omega \right)}=\frac{{\theta }_0{\left(\mathrm{i}\omega \right)}^q+{\theta }_1{\left(\mathrm{i}\omega \right)}^{q-1}+\cdots +{\theta }_q}{{\varphi }_0{\left(\mathrm{i}\omega \right)}^p+{\varphi }_1{\left(\mathrm{i}\omega \right)}^{p-1}+\cdots +{\varphi }_p},$$

(22)

which is described as the frequency response function of the system, arises when the forcing function $\zeta \left(t\right)$ of Equation (18) is replaced by the complex exponential function $exp\left\{\mathrm{i}\omega t\right\}$. It is just the transfer function evaluated by substituting $s=\mathrm{i}\omega .$

In describing the mapping from the forcing function $\zeta \left(t\right)$ to the output variable $y\left(t\right)$, we may consider the rational form of the differential equation, which, on the assumption that there are no repeated roots, has the following partial-fraction decomposition:

$$\begin{array}{ccc}y\left(t\right)=\frac{\theta \left(D\right)}{\varphi \left(D\right)}\zeta \left(t\right)\hfill & =& \left\{\frac{c_1}{D-{\kappa }_1}+\frac{c_2}{D-{\kappa }_2}+\cdots +\frac{c_p}{D-{\kappa }_p}\right\}\zeta \left(t\right)\hfill \\ & =& {\int }_0^{\infty }\left\{c_1{e}^{{\kappa }_1\tau }+c_2{e}^{{\kappa }_2\tau }+\cdots +c_p{e}^{{\kappa }_p\tau }\right\}\zeta (t-\tau )d\tau .\hfill \end{array}$$

(23)

The stability of the system requires all of the values ${\kappa }_j$ to lie in the left half of the complex plane.

The sequence $\{y_t;t=0\pm 1,\pm 2,\dots \}$, which has been sampled at unit intervals for the continuous trajectory of $y\left(t\right)$, can be construed as a process generated by a discrete-time ARMA$(p,p-1)$ model. The latter is the discrete counterpart of the LSDE$(p,q)$ model, and it is described as the exact discrete linear model (EDLM). This result has been demonstrated by Chambers and Thornton (2012) within the context of a state-space representation.

The EDLM can be derived, in theory, by converting each partial-fraction component of the LSDE into a corresponding discrete-time component. In the case of the jth component, which may be real or complex valued, there is

$${\upsilon }_j\left(t\right)=\frac{1}{D-{\kappa }_j}\zeta \left(t\right)={\int }_0^{\infty }{e}^{{\kappa }_j\tau }\zeta (t-\tau )d\tau ={\int }_{-\infty }^t{e}^{{\kappa }_j(t-\tau )}\zeta \left(\tau \right)d\tau .$$

(24)

This result is to be found in Priestley (1981, p. 164) and elsewhere. Within the final expression, the integral on the interval $(-\infty ,t]$ can be separated into two parts, which are the integrals over $(-\infty ,t-1]$ and $(t-1,t]$:

$$\begin{array}{ccc}\hfill {\upsilon }_j\left(t\right)& =& e^{{\kappa }_j}{\int }_{-\infty }^{t-1}{e}^{{\kappa }_j(t-1-\tau )}\zeta \left(\tau \right)d\tau +{\int }_{t-1}^t{e}^{{\kappa }_j(t-\tau )}\zeta \left(\tau \right)d\tau \hfill \\ & =& {\mu }_j{\upsilon }_j(t-1)+{\epsilon }_j\left(t\right),\hfill \end{array}$$

(25)

where ${\mu }_j=e^{{\kappa }_j}$, and where ${\epsilon }_j\left(t\right)$ is a discrete white-noise process. This can be written as $(1-{\mu }_{j}L){\nu }_j\left(t\right)={\epsilon }_j\left(t\right)$, where L is the lag operator that has the effect that $Lx\left(t\right)=x(t-1)$ when applied to a sequence $x\left(t\right)=\{x_t;t=0,\pm 1,\pm 2,\dots \}$.

The assemblage of these discrete-time process ${\upsilon }_j\left(t\right);j=1,2,\dots ,p$ gives rise to an ARMA$(p,p-1)$ model in the form

$$\begin{array}{ccc}\hfill y\left(t\right)=\frac{\beta \left(L\right)}{\alpha \left(L\right)}\epsilon \left(t\right)& =& \left\{\frac{d_1}{1-{\mu }_{1}L}+\frac{d_2}{1-{\mu }_{2}L}+\cdots +\frac{d_p}{1-{\mu }_{p}L}\right\}\epsilon \left(t\right)\hfill \\ & =& \sum _{\tau =0}^{\infty }\left\{d_1{\mu }_1^{\tau }+d_2{\mu }_2^{\tau }+\cdots +d_p{\mu }_p^{\tau }\right\}\epsilon (t-\tau ),\hfill \end{array}$$

(26)

where

$$\begin{array}{ccc}\alpha \left(L\right)\hfill & =& {\alpha }_0+{\alpha }_{1}L+\cdots +{\alpha }_p{L}^p=\prod _{i=1}^p(1-{\mu }_{i}L)\mathrm{and}\hfill \\ \beta \left(L\right)\hfill & =& {\beta }_0+{\beta }_{1}L+\cdots +{\beta }_q{L}^q=\sum _{j=1}^p{d}_j\prod _{i\ne j}^p(1-{\mu }_{i}L).\hfill \end{array}$$

(27)

The white-noise sequence $\epsilon \left(t\right)$ is formed from a weighted combination of the white-noise sequences ${\epsilon }_j\left(t\right);j=1,\dots ,p$. Equation (26) represents the EDLM that is the discrete-time counterpart of the LSDE.

A straightforward comparison of the continuous-time LSDE and the discrete-time EDLM is via their impulse response functions.

The function

$$\psi \left(\tau \right)=c_1{e}^{{\kappa }_1\tau }+c_2{e}^{{\kappa }_2\tau }+\cdots +c_p{e}^{{\kappa }_p\tau },$$

(28)

which is found within Equation (23), is the impulse response function that characterises the mapping of the continuous LSDE model from the forcing function to the output. In the engineering literature, such functions are commonly known as Green’s functions.

According to the sifting property of Dirac’s delta function $\delta \left(\tau \right)$, there is

$$\psi \left(t\right)={\int }_{-\infty }^{\infty }\psi (t-\tau )\delta \left(\tau \right)d\tau ={\int }_{-\infty }^{\infty }\psi \left(\tau \right)\delta (t-\tau )d\tau .$$

(29)

Thus, the impulse response function can be obtained by replacing the forcing function $\zeta \left(t\right)$ in (18) by $\delta \left(t\right)$ and by extending the integral over the negative real line, if necessary, depending on the construction of the limit process that leads to $\delta \left(t\right)$.

The analogous impulse response function of the discrete-time model with $p>q$, which represents its response to a discrete-time unit impulse, is

$$\varphi \left(\tau \right)=d_1{\mu }_1^{\tau }+d_2{\mu }_2^{\tau }+\cdots +d_p{\mu }_p^{\tau },$$

(30)

found within Equation (26), for which the condition of convergence is that all of the values ${\mu }_j$ must lie within the unit circle of the complex plane.

A special case of the continuous-time model arises when $\zeta \left(t\right)$ is replaced by a sequence of Dirac deltas separated by unit intervals and weighted by the ordinates of a discrete-time white-noise process. Then, of course, the continuous-time process coincides with the discrete-time process and, within Equation (26), there is $d_j=c_j$ for all j.

An equivalence at the integer points of the impulse response function of an ARMA model and that of a linear stochastic differential equation driven by a frequency-limited white noise will be used in finding a frequency-limited continuous-time version of the ARMA process. The later will be described as a frequency-limited CARMA process, notwithstanding the common use of this acronym to denote an LSDE that is driven by a white-noise process of unbounded frequencies. In this paper, the acronym LSDE will refer, in most instances, to an equation driven by the increments of a Wiener process.

4. Autocovariance Functions and Spectra

In finding the parameters of a unbounded LSDE when its exact discrete counterpart is available, we shall make use of the autocovariance functions of the discrete and continuous processes.

The autocovariance generating function of the ARMA process is given by

$${\gamma }^d\left(z\right)={\sigma }_{\epsilon }^2\frac{\beta \left(z\right)\beta \left(z^{-1}\right)}{\alpha \left(z\right)\alpha \left(z^{-1}\right)}={\sigma }_{\epsilon }^2\varphi \left(z\right)\varphi \left(z^{-1}\right),$$

(31)

where $\alpha \left(z\right)={\alpha }_0+{\alpha }_{1}z+\cdots +{\alpha }_p{z}^p$ and $\beta \left(z\right)={\beta }_0+{\beta }_{1}z+\cdots +{\beta }_q{z}^q$ are the z-transforms of the autoregressive and moving-average operators, respectively, where z is a complex number and where $\varphi \left(z\right)$, which represents the series expansion of the rational function of $\beta \left(z\right)/\alpha \left(z\right)$, is the z-transform of the impulse response function. This may be described as the discrete-time transfer function, and it becomes the frequency response function when $z=exp\left\{\mathrm{i}\omega \right\}$.

Setting $z=exp\left\{\mathrm{i}\omega \right\}$ places it on the unit circle in the complex plane. By running the argument $\omega$ over the interval $(-\pi ,\pi ]$, the spectral density function or “spectrum” $f\left(\omega \right)$ is generated:

$$f\left(\omega \right)=\frac{1}{2\pi }\sum _{\tau =-\infty }^{\infty }{\gamma }_{\tau }{e}^{-\mathrm{i}\omega \tau }=\frac{1}{2\pi }\left\{{\gamma }_0+2\sum _{\tau =1}^{\infty }{\gamma }_{\tau }cos\left(\omega \tau \right)\right\}.$$

(32)

The spectrum is just the discrete-time Fourier transform of the sequence of autocovariances.

The spectrum and the autocovariance function provide alternative and equivalent characterisations of the ARMA process. The autocovariances are recovered from the spectrum via the inverse transform:

$${\gamma }_{\tau }={\int }_{-\pi }^{\pi }{e}^{i\omega \tau }f\left(\omega \right)d\omega .$$

(33)

On defining $\delta \left(z^{-1}\right)={\sigma }_{\epsilon }^2\beta \left(z^{-1}\right)/\alpha \left(z^{-1}\right)$, Equation (31) can be written as

$${\gamma }^d\left(z\right)\alpha \left(z\right)=\delta \left(z^{-1}\right)\beta \left(z\right),$$

(34)

wherein the coefficients of $\delta \left(z^{-1}\right)$ can be found by solving the equations

$$\delta \left(z\right)\alpha \left(z\right)={\sigma }_{\epsilon }^2\beta \left(z\right).$$

(35)

Then, given the known values on the RHS of (31), those equations can be solved for ${\gamma }_0,\dots ,{\gamma }_p$. Thereafter, the equation

$${\alpha }_0{\gamma }_{\tau }+{\alpha }_1{\gamma }_{\tau -1}+\cdots +{\alpha }_p{\gamma }_{\tau -p}=0;\tau >p,$$

(36)

can be solved recursively to provide the succeeding values $\{{\gamma }_{p+1},{\gamma }_{p+2},\dots \}$.

An analytic expression for the autocovariance function of an ARMA$(p,q)$ model with $p>q$ is also available that exploits the expression of (30) for the discrete-time impulse response function:

$$\begin{array}{ccc}\hfill {\gamma }^d\left(\tau \right)& =& {\sigma }_{\epsilon }^2\sum _{j=0}^{\infty }\left(\sum _{k=1}^p{d}_k{\mu }_k^j\right)\left(\sum _{k=1}^p{d}_k{\mu }_k^{j+\tau }\right)\hfill \\ & =& {\sigma }_{\epsilon }^2\sum _{i=1}^p\left\{\sum _{j=1}^p\frac{d_i{d}_j}{1-{\mu }_i{\mu }_j}\right\}{\mu }_i^{\tau }.\hfill \end{array}$$

(37)

The autocovariance function of the continuous-time LSDE process is also found via its impulse response function. It is assumed that

$$E\left\{\zeta \left(t\right)\zeta (t-\tau )\right\}=\delta \left(\tau \right){\sigma }_{\zeta }^2.$$

(38)

Then,

$$\begin{array}{ccc}\hfill {\gamma }^c\left(\tau \right)& =& E\left\{y\right(t\left)y\right(t-\tau \left)\right\}\hfill \\ & =& E\left\{{\int }_0^{\infty }\psi \left(u\right)\zeta (t-u)du{\int }_0^{\infty }\psi \left(v\right)\zeta (t-\tau -v)dv\right\}\hfill \\ & =& {\sigma }_{\zeta }^2{\int }_0^{\infty }\psi \left(v\right)\psi (v+\tau )dv.\hfill \end{array}$$

(39)

Substituting the expression of (28) for the continuous-time impulse response function $\psi \left(t\right)$ into Equation (39) gives

$$\begin{array}{ccc}\hfill {\gamma }^c\left(\tau \right)={\sigma }_{\zeta }^2{\int }_0^{\infty }\psi \left(t\right)\psi (t+\tau )dt& =& {\sigma }_{\zeta }^2\sum _i\sum _j\left\{c_i{c}_j{\int }_0^{\infty }{e}^{({\kappa }_i+{\kappa }_j)t+{\kappa }_i\tau }dt\right\}\hfill \\ & =& {\sigma }_{\zeta }^2\sum _i\left\{\sum _j{c}_i{c}_j\frac{-e^{{\kappa }_i\tau }}{{\kappa }_i+{\kappa }_j}\right\}.\hfill \end{array}$$

(40)

This expression, which is liable to contain complex-valued terms, may be rendered in real terms by coupling the various conjugate complex terms.

The spectral density function $f\left(\omega \right)$ of the continuous-time process is the Fourier integral transform of autocovariance function:

$$f\left(\omega \right)=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }{\gamma }^c\left(\tau \right)e^{-\mathrm{i}\omega \tau }d\tau =\frac{1}{\pi }{\int }_0^{\infty }{\gamma }^c\left(\tau \right)cos\left(\omega \tau \right)d\tau .$$

(41)

This is a symmetric function; but, in contrast to the autocovariance function of an ARMA process, it is not a periodic function. In the case of an LSDE$(p,q)$ model of Equation (18), the spectral density function is given by

$$f\left(\omega \right)=\frac{{\sigma }_{\zeta }^2}{2\pi }\frac{|{\theta }_0{\left(\mathrm{i}\omega \right)}^p+{\theta }_1{\left(\mathrm{i}\omega \right)}^{p-1}+\cdots +{\theta }_q{|}^2}{{|\left(\mathrm{i}\omega \right)}^p+{\varphi }_1{\left(\mathrm{i}\omega \right)}^{p-1}+\cdots +{\varphi }_p{|}^2}.$$

(42)

This embodies the squared modulus of the complex frequency response function of (22). The connection between the formulations of (41) and (42) is by virtue of the Wiener–Khintchine theorem concerning the spectral representation of a covariance stationary process. This result is to be found in Priestley (1981, p. 250) and elsewhere.

Example 1.

The rational transfer function of the LSDE(2, 1) model with complex-valued poles is represented by

$$\frac{{\theta }_{0}s+{\theta }_1}{s^2+{\varphi }_{1}s+{\varphi }_2}=\frac{c}{s-\kappa }+\frac{c^{*}}{s-{\kappa }^{*}}=\frac{c(s-{\kappa }^{*})+c^{*}(s-\kappa )}{(s-\kappa )(s-{\kappa }^{*})},$$

(43)

where $\kappa =\delta +\mathrm{i}\omega$, ${\kappa }^{*}=\delta -\mathrm{i}\omega$ are the poles and where $c=a+\mathrm{i}b$, $c^{*}=a-\mathrm{i}b$ are the numerator coefficients of the partial fraction expansion. It follows that

$$\begin{array}{ccc}\hfill {\theta }_0& =& 2aand{\theta }_1=2(a\delta +b\omega )and,conversely,\hfill \\ \hfill a& =& \frac{{\varphi }_0}{2}andb=\frac{{\varphi }_0\delta +{\theta }_1}{2}.\hfill \end{array}$$

(44)

According to the formula of (40), the autocovariance function of the LSDE(2, 1) process, is given by

$$\begin{array}{ccc}\hfill {\gamma }^c\left(\tau \right)& =& -{\sigma }_{\zeta }^2\left\{\frac{c^2}{2\kappa }+\frac{c{c}^{*}}{\kappa +{\kappa }^{*}}\right\}{e}^{\kappa \tau }-{\sigma }_{\zeta }^2\left\{\frac{{\left(c^{*}\right)}^2}{2{\kappa }^{*}}+\frac{c^{*}c}{{\kappa }^{*}+\kappa }\right\}{e}^{{\kappa }^{*}\tau }\hfill \\ & =& -{\sigma }_{\zeta }^2\left\{\frac{c^2{\kappa }^{*}{e}^{\kappa \tau }+{\left(c^{*}\right)}^2\kappa e^{{\kappa }^{*}\tau }}{2\kappa {\kappa }^{*}}\right\}-{\sigma }_{\zeta }^2\frac{c{c}^{*}}{k+k^{*}}\left(e^{\kappa \tau }+e^{{\kappa }^{*}\tau }\right).\hfill \end{array}$$

(45)

Here, the second expression gathers the conjugate complex terms. The conjugate complex exponential functions will give rise to real-valued trigonometric functions, which will embody the following terms:

$$\begin{array}{ccc}\hfill \kappa {\kappa }^{*}& =& {\delta }^2+{\omega }^2,\kappa +{\kappa }^{*}=2\delta c{c}^{*}=a^2+b^2,\hfill \\ \hfill c^2{\kappa }^{*}& =& {(a+\mathrm{i}b)}^2(\delta -\mathrm{i}\omega )=\delta (a^2-b^2)+2ab\omega +\mathrm{i}\left\{2ab\delta -\omega ({\alpha }^2+b^2)\right\},\hfill \\ \hfill {\left(c^{*}\right)}^2\kappa & =& {(a-\mathrm{i}b)}^2(\delta +\mathrm{i}\omega )=\delta (a^2-b^2)+2ab\omega +\mathrm{i}\left\{\omega (a^2-b^2)-2ab\delta \right\}.\hfill \end{array}$$

(46)

The terms may be combined using the result that, if $\kappa =\delta +\mathrm{i}\omega$, ${\kappa }^{*}=\delta -\mathrm{i}\omega$ and $c=a+\mathrm{i}b$, $c^{*}=a-\mathrm{i}b$, then

$$\psi \left(t\right)=c{e}^{\kappa \tau }+c^{*}{e}^{{\kappa }^{*}\tau }=2e^{\delta \tau }\{acos\left(\omega \tau \right)-bsin\left(\omega \tau \right)\}.$$

(47)

It follows that an alternative expression for the autocovariance function is

$$\begin{array}{ccc}\hfill {\gamma }^c\left(\tau \right)=\frac{{\sigma }_{\zeta }^2}{{\delta }^2+{\omega }^2}{e}^{\delta \tau }\left\{\right[2ab\delta & -& \omega (a^2+b^2)]sin\left(\omega \tau \right)-[2ab\omega +\delta (a^2-b^2)]cos\left(\omega \tau \right)\}\hfill \\ & -& {\sigma }_{\zeta }^2{e}^{\delta \tau }\frac{(a^2+b^2)}{\delta }cos\left(\omega \tau \right).\hfill \end{array}$$

(48)

There is no difficulty in dealing with models of higher orders.

The expression of (48) for the autocovariance function of an LSDE(2, 1) process with conjugate complex poles will be employed in a subsequent example.

5. The Continuous-Time Frequency-Limited ARMA Process

The basis of the frequency-limited continuous-time CARMA model is a frequency-limited white-noise forcing function. This is derived from a discrete-time white-noise process $\epsilon \left(t\right)=\{{\epsilon }_t;t=0,\pm 1,\pm 2,\dots \}$, in accordance with the Nyquist–Shannon sampling theorem, by the simple expedient of attaching a sinc function to each of its ordinates and by adding together the overlapping kernel functions. The resulting process may be denoted by

$${\epsilon }_c\left(t\right)=\sum _{k=-\infty }^{\infty }{\epsilon }_k\phi (t-k),\mathrm{with}t\in \mathcal{R},$$

(49)

where $\phi (t-k)=sin\left\{\pi \right(t-k\left)\right\}/\pi (t-k)$ is the sinc function centred on $t=k$. Hereafter, we shall omit the subscript c, which serves to distinguish the continuous function ${\epsilon }_c\left(t\right);t\in \mathcal{R}$ from the sequence $\epsilon \left(t\right);t\in \{0,\pm 1,\pm 2,\dots \}=\mathcal{Z}$.

The continuous frequency-limited white noise, which is assumed to have a mean of zero, has an autocovariance function that is a scaled sinc function:

$${\gamma }_{\epsilon }\left(\tau \right)={\sigma }_{\epsilon }^2\phi \left(\tau \right).$$

(50)

To understand this, consider two elements of the process. These are ${\epsilon }_t$, which is the amplitude of the sinc function located at t, and ${\epsilon }_s={\epsilon }_t\phi (t-s)+\upsilon$. The latter consists of the attenuated tail of the sinc function at t and a term, $\upsilon$, representing the combined effects of the sinc functions at all other locations. Since $\upsilon$ is independent of ${\epsilon }_t$ it follows that

$${\gamma }_{\epsilon }\left(\tau \right)=E\left({\epsilon }_t{\epsilon }_s\right)={\sigma }_{\epsilon }^2\phi (t-s);\tau =t-s.$$

(51)

It may be observed that $\phi \left(z\right)$ is symmetric, with $\phi \left(z\right)=\phi \left(z^{-1}\right)$ and idempotent, with $\phi \left(z\right)={\phi }^2\left(z\right)$. Therefore, $\phi \left(z\right)=\phi \left(z^{-1}\right)\phi \left(z\right)$ is its own autocovariance generating function. The condition that $\phi \left(\tau \right)=0$ for $\tau \in \{\pm 1,\pm 2,\dots \}$ indicates that the set of sinc functions displaced one from another by multiples of one sample period form an othogonal basis. This corresponds to the fact that a continuous white-noise process with a frequency limit of $\pi$ will deliver a discrete-time white-noise sequence when sampled at unit intervals. When sampled at other intervals, the resulting sequence will manifest some autocorrelation.

On defining $y\left(t\right)={\sum }_k{y}_k\psi (t-k)$, the equation of the frequency-limited continuous-time ARMA model may be denoted by

$$\sum _{j=0}^p{\alpha }_{j}y(t-j)=\sum _{j=0}^q{\beta }_j\epsilon (t-j),$$

(52)

with ${\alpha }_0=1$. The model has a moving-average representation in the form of

$$y\left(t\right)=\sum _{j=0}^{\infty }{\varphi }_j\epsilon (t-j),$$

(53)

where the coefficients are from the series expansion of the rational function $\beta \left(z\right)/\alpha \left(z\right)=\varphi \left(z\right)$. In this notation, the frequency-limited continuous-time ARMA model is indistinguishable from the discrete-time model.

If $y\left(t\right)={\sum }_i{\varphi }_i\epsilon (t-i)$, and $y\left(s\right)={\sum }_j{\varphi }_j\epsilon (s-j)$, with $t,s\in \mathcal{R}$ and $i,j\in \mathcal{Z}$, then the autocovariance function of $y_t=y\left(t\right)$ and $y_s=y\left(s\right)$ is

$$\begin{array}{ccc}\hfill C\left\{\sum _i{\varphi }_i{\epsilon }_{t-i},\sum _j{\varphi }_j{\epsilon }_{s-j}\right\}& =& \sum _i\sum _k{\varphi }_i{\varphi }_{i+k}C({\epsilon }_t,{\epsilon }_{s-k});k=j-i\hfill \\ & =& \sum _k{\gamma }_k\phi (\tau -k)=\gamma \left(\tau \right);\tau =t-s,\hfill \end{array}$$

(54)

where C stands for the covariance operator and where ${\gamma }_k={\sigma }_{\epsilon }^2{\sum }_j{\varphi }_j{\varphi }_{j+j}$ is the kth autocovariance of the discrete-time process. It can be seen immediately that $\gamma \left(\tau \right)={\gamma }_{\tau }$, is a discrete-time autocovariance when $\tau$ takes an integer value, and that the continuous-time autocovariance function is obtained from the discrete-time function by sinc-function interpolation.

A difficulty in realising this expression for the autocovariance function lies in the fact that the sinc functions are supported on the entire real line. An alternative approach is to pursue the (inverse) Fourier transform of the spectral density function. This function is equivalent to the central segment of the periodic spectral density function $f\left(\omega \right)$ of the discrete-time ARMA model.

Thus,

$$\gamma \left(\tau \right)={\int }_{-\pi }^{\pi }{e}^{\mathrm{i}\omega \tau }f\left(\omega \right)d\omega ={\int }_0^{\pi }2cos\left(\omega \tau \right)f\left(\omega \right)d\omega ,$$

(55)

where the second equality follows in consequence of the symmetry of $f\left(\omega \right)=f(-\omega )$. To realise this expression it is necessary to approximate the Fourier integral by a discrete cosine Fourier transform embodying a large the number of points sampled from the function $f\left(\omega \right)$ over the interval $[-\pi ,\pi ]$.

Although such an approximation of the integral can be highly accurate, it is altogether easier to generate the continuous autocovariance function by allowing the argument $\tau$ in the analytic expression of (37) for the discrete-time autocovariance function to vary continuously.

6. Estimates of the Linear Stochastic Models

The available computer programs for estimating ARMA models are numerous and well developed. Those for estimating LSDE models are comparatively rare, and some of them are restricted to continuous-time autoregressive models that lack moving-average components. All of them presuppose that the processes are driven by forcing functions of unlimited frequencies.

The direct approach to the estimation of a linear stochastic differential equation (LSDE) depends on the optimisation of a criterion function such as the likelihood function or a residual sum of squares. Typically, a state-space representation of the model is employed in evaluating the function and in seeking its optimum.

The early approaches to estimating an LSDE concentrated on pure autoregressive models that lack a moving-average component. Thus, Bergstrom (1966) adopted purely autoregressive formulations in constructing multi-equation continuous-time econometric models. An interesting and a unique approach to the estimation of a single purely autoregressive equation has been taken, more recently, by Hyndman (1993).

An important contribution to the estimation via the direct approach was made by Jones (1981), who provided a computer program for estimating an autoregressive LSDE that was subject to a white-noise contamination. More recent developments, which allow for a moving average component and for the presence of stochastic trends, have been due to Harvey and Stock (1988, 1985) and to Chambers and Thornton (2012).

In the alternative indirect approach, the LSDE estimates are derived from prior estimates of an autoregressive moving-average (ARMA) model. First, the autoregressive dynamics of the ARMA model are translated to the continuous-time dynamics of the LSDE via a one-to-one mapping. Then, a moving-average component is sought that will match the continuous-time autocovariance function to the discrete-time function at the integer lags. An ARMA$(p,q)$ model with $q<p$ will lead, invariably, to an LSDE$(p,p-1)$ model.

The indirect approach has received less attention than the direct approach. However, Söderström (1991) has surveyed some of the available methods and he has developed a state-space method that has been appraised, in comparison with other methods, by Larsson et al. (2006).

The methods of this paper also presuppose the availability of valid ARMA estimates to be used as the starting point. In the case of an ARMA model that is deemed to be free of the effects of aliasing, the corresponding frequency-limited continuous-time CARMA model is available through a one-to-one correspondence based on the equivalence of the impulse response functions.

In some cases, the estimates of an ARMA$(p,p-1)$ model are to be regarded as those of an exact discrete linear model, or EDLM, that is the counterpart of a linear stochastic differential equation, or LSDE, driven by a forcing function of unbounded frequencies.

In such cases, the EDLM will suffer the effects of aliasing, and the continuous-time moving-average parameters are no longer readily available from their discrete-time counterparts. Then, a special iterative estimation procedure is called for that fulfils the principle of autocovariance equivalence.

The task may also arise of converting an LSDE driven by a forcing function of unbounded frequencies to its exact discrete-time counterpart model, or EDLM. This is bound to be an ARMA$(p,p-1)$ model, in the case of an LSDE$(p,q)$ model with $p>q$.

The principle of autocovariance equivalence indicates that this conversion may be achieved by applying a Cramér–Wold decomposition to the autocovariance function of the LSDE, which is shared by the EDLM. Since the discrete-time autoregressive parameters can be inferred directly from their continuous-time counterparts, the Cramér–Wold decomposition is concerned only with finding the moving-average parameters.

6.1. From ARMA to CARMA

The frequency-limited CARMA model is the continuous-time version of the discrete-time ARMA model. To find a continuous trajectory that interpolates the data from which the ARMA estimates have been derived, it is sufficient to replace the discrete temporal index of Equation (2) by a continuous variable $t\in [0,T)$. This amounts to a Fourier interpolation, which is equivalent to a sinc function interpolation employing Dirichelet kernels.

To find a differential equation representing the frequency-limited CARMA model, it is necessary to find the partial-fraction decomposition of the ARMA operator $\beta \left(z\right)/\alpha \left(z\right)$:

$$\frac{{\beta }_0+{\beta }_{1}z+\cdots +{\beta }_q{z}^q}{(1-{\mu }_{1}z)(1-{\mu }_{2}z)\cdots (1-{\mu }_{p}z)}=\frac{d_1}{1-{\mu }_{1}z}+\frac{d_2}{1-{\mu }_{2}z}+\cdots +\frac{d_p}{1-{\mu }_{p}z}.$$

(56)

The jth partial fraction coefficient is

$$d_j=\frac{{\beta }_0+{\beta }_1/{\mu }_j+\cdots +{\beta }_q/{\mu }_j^q}{{\prod }_{i\ne j}(1-{\mu }_i/{\mu }_j)}={\mu }_j^{p-1-q}\frac{({\beta }_0{\mu }_j^q+{\beta }_1{\mu }_j^{q-1}+\cdots +{\beta }_q)}{{\prod }_{i\ne j}({\mu }_j-{\mu }_i)}.$$

(57)

The roots ${\mu }_1,{\mu }_2,\dots ,{\mu }_p$ of the polynomial $\alpha \left(z\right)$ can be calculated reliably via the procedure of Müller (1956), of which versions has been coded by Pollock (1999) in Pascal and C. The numerator coefficients $d_1,d_2,\dots ,d_p$ are best calculated via the expression on the RHS of (57).

A complex pole of the ARMA model (i.e., a root of the equation $\alpha \left(z\right)=0$) takes the form of

$$\mu =\alpha +\mathrm{i}\beta =\rho \{cos\left(\omega \right)+\mathrm{i}sin\left(\omega \right)\}=\rho e^{\mathrm{i}\omega },$$

(58)

where

$$\rho =\sqrt{{\alpha }^2+{\beta }^2}\mathrm{and}\omega ={tan}^{-1}\left(\frac{\beta }{\alpha }\right),$$

(59)

Here, it is assumed that $\omega \in [0,\pi ]$ and that $\rho \in (0,1)$. The corresponding pole of the frequency-limited CARMA differential equation is

$$\kappa =ln\left(\mu \right)=ln\left(\rho \right)+\mathrm{i}\omega =\delta +\mathrm{i}\omega ,$$

(60)

with $\delta \in (-\infty ,0)$, which puts it in the left half of the s-plane, as is necessary for the stability of the system.

The remaining task is to assemble the roots and the coefficients to form the rational operator $\theta \left(s\right)/\varphi \left(s\right)$ of the differential equation:

$$\frac{{\theta }_0{s}^q+{\theta }_1{s}^{q-1}+\cdots +{\theta }_q}{{\varphi }_0{s}^p+{\varphi }_1{s}^{p-1}+\cdots +{\varphi }_p}=\frac{d_1}{s-{\kappa }_1}+\frac{d_2}{s-{\kappa }_2}+\cdots +\frac{d_p}{s-{\kappa }_p}.$$

(61)

In effect, the differential equation of the frequency-limited CARMA model is provided by Equation (23), with $c_j=d_j$ for $j=1,2,\dots ,p$.

Given the latter equalities and the correspondence between the ARMA poles ${\mu }_j$ and the poles ${\kappa }_j=ln\left({\mu }_j\right)$ of the differential equation, it follows that the impulse response function of the frequency-limited CARMA model, provided by Equation (28), is equal, at the integer points, to the impulse response function of the ARMA model, provided by Equation (30). Likewise, the continuous autocovariance function of the CARMA model with be equal, at the integer points, to that of the ARMA model. Therefore, in this case, impulse invariance is equivalent to autocovariance equivalence.

The principal of impulse invariance is familiar in communications engineering as a means of deriving a discrete-time filter from its continuous-time counterpart. See, for example, Oppenheim and Schafer (2013).

Example 2.

To illustrate the mapping from the discrete-time ARMA model to a continuous-time frequency-limited CARMA model, an ARMA(2, 1) model is chosen with conjugate complex poles $\alpha \pm \mathrm{i}\beta =\rho exp\{\pm \mathrm{i}\theta \}$, where $\rho =\sqrt{{\alpha }^2+{\beta }^2}=0.9$ and $\theta ={tan}^{-1}(\beta /\alpha )=\pi /4={45}^{\circ }$. The moving-average component has a zero of $0.5$. The ARMA process generates prominent cycles of an average duration of roughly 8 periods.

The spectral density function of the ARMA process is illustrated in Figure 1. Here, it will be observed that the function is virtually zero at the limiting Nyquist frequency of π. Therefore, it is reasonable to propose that the corresponding continuous-time model should be driven by a white-noise forcing function that is bounded by the Nyquist frequency.

The parameters of the resulting continuous-time frequency-limited CARMA model are displayed below, beside those of the ARMA model:

$$\begin{array}{c}ARMA\hfill \\ {\alpha }_0=1.0\hfill \\ {\alpha }_1=-1.2728\hfill \\ {\alpha }_2=0.8100\hfill \\ \hfill \\ {\beta }_0=1.0\hfill \\ {\beta }_1=-0.5\hfill \end{array}\begin{array}{c}CARMA\hfill \\ {\varphi }_0=1.0\hfill \\ {\varphi }_1=0.2107\hfill \\ {\varphi }_2=0.6280\hfill \\ \hfill \\ {\theta }_0=1.0\hfill \\ {\theta }_1=0.2737\hfill \end{array}$$

In Figure 2, the discrete autocovariance function of the ARMA process is superimposed on the continuous autocovariance function of the frequency-limited CARMA process. The former has been generated by the procedure described by the recursive Equations (34)–(36). The latter has been generated by Equation (37), wherein the index τ varies continuously.

The spectral density function of the frequency-limited CARMA process is the integral Fourier transform of the continuous autocovariance function, whereas the spectral density function of the ARMA process is the discrete Fourier transform of the autocovariance sequence. The frequency limitation of the CARMA process means that there is no aliasing in the sampling process. Therefore, the two spectra are identical within the Nyquist interval.

6.2. From ARMA to LSDE

When a linear stochastic differential equation (LSDE) is driven by a white-noise forcing function of unbounded frequencies, the ARMA model that represents its discrete-time counterpart, which is the EDLM, is bound to be affected by aliasing, even if this is to a negligible extent. The aliasing also affects the impulse response function. Therefore, the principle of impulse invariance can no longer serve to determine the parameters of the continuous-time model. Instead, the principle of autocovariance equivalence, enunciated by Bartlett (1946), must be relied upon to translate from the discrete-time model to the continuous-time model.

This principal asserts that the continuous-time autocovariance function must be equal, at the integer points, to that of the discrete-time model. Given that the autoregressive parameters of the LSDE are implied by those of the EDLM, it follows that only the moving-average parameters are liable to be determined in fulfilment of the principle.

Let ${\gamma }^d(\mu ,d)$ be the autocovariance function of the EDLM, as given by Equation (37), and let ${\gamma }^c(\kappa ,c)$ be the corresponding autocovariance function of the LSDE, which is sampled at the integer points from the function of (40). Since Equation (60) expresses the poles of the LSDE in terms of those of the ARMA model, there is $\kappa =\kappa \left(\mu \right)$, and the principle of autocovariance equivalence can be expressed via the equation

$${\gamma }^c\{\kappa \left(\mu \right),c\}={\gamma }^d(\mu ,d).$$

(62)

Then, the parameters of the LSDE can be derived once a value of $c=[c_1,c_2,\dots ,c_p]$ has been found that satisfies this equation.

The value of c can be found by using an iterative optimisation procedure to find the zero of the function

$$z\left(c\right)=\sum _{\tau =0}^p{w}_{\tau }{\{{\gamma }_{\tau }^c\left(c\right)-{\gamma }_{\tau }^d\}}^2.$$

(63)

wherein $w_1,w_2,\dots ,w_p$ are a set of positive weights. The weights should normally be set to unity, but the use of differential weights might assist the convergence of the iterations.

It will be observed that the coefficients $c_j$ of the continuous autocovariance function of (40) are liable to be complex-valued. This can make it tedious, if not difficult, to calculate the analytic derivatives of $z\left(c\right)$, which might be required by the optimisation procedure. Therefore, the derivative-free procedure of Nelder and Mead (1965) can be used to advantage. The algorithm has been coded in Pascal by Bundy and Garside (1978).

The autocovariance function of the LSDE can also be rendered in a form suggested by Equation (48), wherein the parameters $c=(a,b)$ are real-valued.

Söderström (1990) has determined the possible locations of the zeros of an EDLM(2, 1) as a function of the poles of the LSDE to which it corresponds; and he has shown that the zeros can only occur in a restricted set within the unit circle. This affects the possibility of finding an LSDE to correspond to a specified ARMA model.

An ingredient of the algorithm of Söderström (1991) for translating an ARMA model to an LSDE is a spectral factorisation. If the factorisation cannot be accomplished, then there is no corresponding continuous-time model. Therefore, it has been proposed that the algorithm is an appropriate means of determining whether there exists an LSDE corresponding to a specified ARMA model.

Whenever there exists a continuous-time LSDE of which the sampled autocovariance function conforms to that of an ARMA process, it is said that the ARMA process is embedded in the LSDE. Brockwell (1995) has also investigated the conditions that must be fulfilled by the ARMA process if it is to be embeddable. He has shown that every discrete-time AR(2) process with distinct complex conjugate roots can be embedded in either a continuous-time LSDE(2, 1) process or a in continuous-time LSDE(4, 2) process, or, in some cases, in both.

For the algorithm that we are pursuing, the non-availability of an LSDE corresponding to the ARMA model is indicted by a nonzero minimum of the criterion function. However, when the minimum is close to zero, then the ARMA model will be close to the EDML that corresponds to the minimising LSDE.

In such cases, we might surmise that the divergence of the ARMA model from the EDLM may be due to the vagaries of the processes of sampling and estimation that have given rise to the ARMA model. We might also surmise that the minimising LSDE will be close to one that would be delivered by a direct method of estimation applied to the sample from which the ARMA model has arisen.

Example 3.

The mapping from the discrete-time ARMA model to a continuous-time LSDE model can be illustrated, in the first instance, with the ARMA(2, 1) model of Example 2.

The parameters of the corresponding LSDE(2, 1) model are obtained by using the Nelder–Mead procedure to find the minimum of the criterion function of (63), where it is assumed that the variance of the forcing function is ${\sigma }_{\zeta }^2=1$.

There are four points that correspond to zero-valued minima, where the ordinates of the discrete and continuous autocovariance functions coincide at the integer points. These four points, together with the corresponding moving-average parameters, are shown in

Table 1. The various values to the moving-average parameters of the LSDE(2, 1) model inferred from the discrete-time autocovariance function according to the principle of autocovariance equivalence. The parameters a and b are arguments of the continuous autocovariance function of (48).

	a	b	${\mathit{\theta }}_0$	${\mathit{\theta }}_1$
(i)	−0.4544	0.2956	−0.9088	0.5601
(ii)	0.4544	0.4175	0.9088	0.5601
(iii)	−0.4544	−0.4174	−0.9088	−0.5601
(iv)	0.4544	−0.2956	0.9088	−0.5601

.

Here, the parameter values of (i) and (iv) are equivalent, as are those of (ii) and (iii). Their difference is a change of sign, which can be eliminated by normalising ${\theta }_0$ at unity and by adjusting variance of the forcing function accordingly.

The miniphase condition, which corresponds to the invertibility condition of a discrete-time model, requires the zeros to be in the left half of the s-plane. Therefore, (ii) and (iii) on the NE–SW axis are the chosen pair.

These estimates of the LSDE(2, 1) are juxtaposed below with those of the frequency-limited CARMA(2, 1) model derived from the same ARMA model:

$$\begin{array}{c}CARMA\hfill \\ {\varphi }_0=1.0\hfill \\ {\varphi }_1=0.2107\hfill \\ {\varphi }_2=0.6280\hfill \\ \hfill \\ {\theta }_0=1.0\hfill \\ {\theta }_1=0.2737\hfill \end{array}\begin{array}{c}LSDE\hfill \\ {\varphi }_0=1.0\hfill \\ {\varphi }_1=0.2107\hfill \\ {\varphi }_2=0.6280\hfill \\ \hfill \\ {\theta }_0=0.9088\hfill \\ {\theta }_1=0.5601\hfill \end{array}$$

The autoregressive parameters of the frequency-limited CARMA model and of the LSDE model are, of course, identical. However, there is a surprising disparity between the two sets of moving-average parameters. Nevertheless, when they are superimposed on the same diagram—which is Figure 3—the spectra of the two models are seen virtually to coincide. Moreover, the parameters of the ARMA model can be recovered exactly from those of the LSDE by an inverse transformation, which will be described later.

The explanation for this outcome is to be found in the remarkable flatness of the criterion function in the vicinity of the minimising points, which are marked on both sides of Figure 4 by black dots. The flatness implies that a wide spectrum of the parameter values of the LSDE will give rise to almost identical autocovariance functions and spectra.

The left side of Figure 4 show the some equally-spaced contours of the z-surface of the criterion function, which are rising from an annulus that contains the minima. The minima resemble small indentations in the broad brim of a hat.

The right side of Figure 4, which is intended to provide more evidence of the nature of the criterion function in the vicinity of the minima, shows the contours of the function $q=1/(z+d)$, where d is a small positive number that prevents a division by zero. We set $d=(X-RM)/(R-1)$, where $M=min\left(z\right)$, $X=max\left(z\right)$ and where $R=max\left(q\right)/min\left(q\right)=60$. The extended lenticular contours surrounding the minima are a testimony to the virtual equivalence of a wide spectrum of parameter values.

A variant to the ARMA(2, 1) model is one that has a pair of complex conjugate poles $\rho exp\{\pm \mathrm{i}\theta \}$ with the same argument as before, which is $\theta ={tan}^{-1}(\beta /\alpha )=\pi /4={45}^{\circ }$, and with a modulus that has been reduced to $\rho =0.5$. The model retains the zero of the transfer function (i.e., the root of the moving-average operator) of $0.5$. The ARMA parameters and those of the corresponding LSDE are as follows:

$$\begin{array}{c}ARMA\hfill \\ {\alpha }_0=1.0\hfill \\ {\alpha }_1=-0.7071\hfill \\ {\alpha }_2=0.2500\hfill \\ \hfill \\ {\beta }_0=1.0\hfill \\ {\beta }_1=-0.5\hfill \end{array}\begin{array}{c}LSDE\hfill \\ {\varphi }_0=1.0\hfill \\ {\varphi }_1=1.3863\hfill \\ {\varphi }_2=1.0973\hfill \\ \hfill \\ {\theta }_0=1.5012\hfill \\ {\theta }_1=0.8905\hfill \end{array}$$

Figure 5 shows the spectral density function of the LSDE and of the ARMA model superimposed on same diagram. The spectrum of the LSDE extends far beyond the Nyquist frequency of π, which is the limiting ARMA frequency.

The ARMA process, which is to be regarded as a sampled version of the LSDE, which is the EDLM, is seen to suffer from a high degree of aliasing, whereby the spectral power of the LSDE that lies beyond the Nyquist frequency is mapped into the Nyquist interval $[-\pi ,\pi ]$, with the effect that the profile of the ARMA spectrum is raised considerably. On this basis, it can be asserted that the ARMA model fails to reveal the frequency content of the underlying continuous-time process. However, its other uses, such as a means of making forecasts, are unaffected.

6.3. From LSDE to ARMA

The mapping from the LSDE to the EDLM also depends on the principle of autocovariance equivalence. Thus, the ARMA parameters will be derived from the sampled ordinates the continuous-time autocovariance function of the LSDE. The essential equation is that of (31), which can be cast in the form of

$$\beta \left(z\right)\beta \left(z^{-1}\right)={\sigma }_{\epsilon }^2\alpha \left(z\right){\gamma }^c\left(z\right)\alpha \left(z^{-1}\right),$$

(64)

where ${\gamma }^c\left(z\right)$ is the z-transform of the elements ${\gamma }_{\tau }^c;\tau \in \{0,\pm 1,\pm 2,\dots \}$ sampled from the autocovariance function to the LSDE.

Since the discrete-time autoregressive parameters within $\alpha \left(z\right)$ can be inferred from those of the LSDE, only the moving-average parameters within $\beta \left(z\right)$ need to be derived from Equation (64). They can be obtained, together with the variance parameter ${\sigma }_{\epsilon }^2$, by a Cramér–Wold Factorisation. A Newton–Raphson algorithm for effecting the factorisation has been presented by Wilson (1969), and it has been coded in Pascal and C by Pollock (1999).

7. Summary and Conclusions

The intention of this paper has been to clarify the relationship between unconditional linear stochastic models in discrete and continuous time, and to provide a secure means of relating such models to each other.

The results of the paper emphasise the importance of an awareness of the frequency-domain characteristics of the forcing functions of such models.

In a recent paper—Pollock (2012)—the author has shown that there are liable to be severe biases in the estimated parameters of ARMA models when the data have a frequency limit significantly below the Nyquist value of $\pi$, which is the maximum frequency for a discrete-time model.

A way of overcoming the problem is to create a continuous data trajectory by Fourier interpolation, from which new values can be resampled at an appropriate rate that is less than the original rate of sampling. The procedure invokes the concept of a frequency-limited continuous-time CARMA model.

The present paper further explores both the frequency-limited CARMA model and the usual linear stochastic differential equations (LSDEs) that are driven by a white-noise forcing functions of unbounded frequencies.

Example 3 shows the aliasing effects that occur when the forcing function has no frequency limit and when the ARMA transfer function imposes only a weak attenuation on the high-frequency elements. This provides a justification for adopting an LSDE of unbounded frequencies for the continuous-time translation of the ARMA model.

Given that continuous-time models are bound to be inferred from discretely sampled data, it is appropriate to base their estimation on the discrete-time autocovariance functions that are provided by ARMA models.

A secure method of translating from a discrete autocovariane function to a continuous model has been proposed. The method has been realised in the computer program CONCRETE.PAS, which is available at the author’s website, where both the compiled program and it code can be found. The program requires the prior specification on an ARMA model in the discrete-time domain, or an LSDE in the continuous-time domain. It proceeds from either of these starting points to find the corresponding model in the other domain. It is intended to develop the program by including direct maximum-likelihood methods of estimation in both the discrete and the continuous domains.

An associated program CONTEXT.PAS is also available, which plots the contour map of the surface of the criterion function that is employed in matching the autocovariance function of the LSDE to that of the ARMA model (Supplementary Materials).

As Söderström (1990, 1991), and others have noted, there are ARMA models for which there are no corresponding LSDE’s. The procedure for translating from an ARMA model to an LSDE that has been described in the present paper reveals such cases by its failure to find a zero-valued minimum of the criterion function. However, it can be relied upon to find the LSDE most closely related to the ARMA model.

One of the bugbears in estimating an LSDE is the need to satisfy the miniphase condition, by ensuring that the zeros fall on the left side of the s-plane. The way to achieve this is to find the zeros by factorising the moving-average polynomial, and then to reflect them across the vertical, imaginary, axis if they happen to fall on the right side. This is an inescapable burden whenever the order of the moving-average component is greater than one.