Construction of an SDE Model from Intraday Copper Futures Prices

Abstract

This paper introduces a model for intraday copper futures prices based on a stochastic differential equation (SDE). In particular, we derive an SDE that fits the model to the data and that is based on the whitening filter approach, a method characterizing linear time-variant systems. This method is applied to construct a model able to simulate the trajectories of copper futures prices, statistically described by means of an empirical autocorrelation approach. We show that the predictability of copper futures prices is rather weak. In fact, the developed model produces trajectories close to the actual data only in the short term. Consequently, the investment risk for copper futures is high. We also show that the performance of the model improves significantly if the time series satisfy particular conditions, e.g., those with a determinism measure.

1. Introduction

Financial time series modeling is an essential aspect of forecasting and risk evaluation in financial markets. In this paper, we consider the case of copper, which, according to the NYMEX, is the third most used metal, which makes it important to assess the nature of its price fluctuations. Further, the likely pressure on copper prices due to its possible scarcity is a cause for concern Gordon et al. (2006); Tilton and Lagos (2007), especially in light of its importance for the growing network industry.

Our previous works, Mastroeni and Vellucci (2019); Mastroeni et al. (2018), showed that a significant level of noise is usually present in the time series of copper futures intraday prices, supporting the conclusion that logarithmic returns have both a stochastic and deterministic nature. Hence, the use of stochastic models appears to be a natural choice.

In this paper, we develop a novel stochastic differential equation (SDE) that models the data and is based on the whitening filter approach Wiener (1949), a method characterizing linear time-variant systems. From a statistical point of view, we assume that the time series of prices can be described by autocorrelation. This property is obtained by means of statistical analyses of the historical data recorded for the futures copper closing prices (HG1 ticker) as exchanged on the COMEX market (CMX). Starting from a model of the empirical autocorrelation shown by the time series of prices, we introduce an SDE, which is characterized by this autocorrelation, and fits the model to the data. To the best of our knowledge, there have been no papers following this approach.

The purpose of the developed SDE model is to move one step further concerning the analysis started in our previous paper Mastroeni et al. (2018). In that paper, we showed that the time series of copper prices (the same studied in the present paper but for a shorter interval) exhibited both stochastic and chaotic features. At the same time, the recurrence plot of the time series revealed a pattern typical of intermittency phenomena.

The presence of stochastic features and intermittency phenomena called for a second step of the time series analysis started in Mastroeni et al. (2018), i.e., a new approach to short-term modeling for copper futures prices.

Our model produced trajectories close to the real data only in the short run. This is due, on the one hand, to the fitted model from the empirical autocorrelation and, on the other, to the predominantly stochastic behavior of the copper price time series, which can only be predicted in the short term Mastroeni and Vellucci (2019, 2022). To investigate the stochastic (vs deterministic) nature of the data, we follow the intuition behind recurrence analysis Eckmann et al. (1987); Marwan et al. (2007). In the Appendix A, we show that the performance of the model improved significantly when the time series had an empirical autocorrelation function close to the damped cosine model and a determinism measure (DET), based on a recurrence plot, close to one. In the numerical section, we show that the presence of a large stochastic component produced a poorer performance than that obtained for time intervals with a large deterministic component.

So far, the literature does not seem to have focused on mathematical modeling of the copper futures market, following econometric and statistical approaches instead. The following is a brief review of the literature supporting this claim.

In Jin et al. (2021), the authors analyzed the effect of the price correlation between domestic and foreign copper futures contracts. Guo et al. (2020) investigated the nonlinear correlation between the spot and futures prices in China’s copper market using nonlinear Granger causality and multifractal methods. Gong and Lin (2018) examined whether structural breaks contained incremental information for forecasting the volatility of copper futures. In Zheng et al. (2022), the authors conducted a comparative exploration of the chaotic characteristics of Chinese and international copper futures prices, showing the differences in those series through recurrence plots and correlation dimensions.

The co-movements of the dynamic correlations between copper futures and spot prices on a scale by scale basis were analyzed in Yu et al. (2021), using grey correlation analysis and wavelet analysis. In Galán-Gutiérrez and Martín-García (2022), Granger causality was exploited to test whether, during the first COVID-19 wave, the evolution of the pandemic was cointegrated with the price structure of copper futures, on the one hand, and with incremental levels of copper stocks, on the other. The authors of Idrovo-Aguirre and Contreras-Reyes (2021) proposed the impulse–response functions of a vector autoregressive model to capture the dynamic between copper prices and house building permits.

As for the development of a model that fit the same copper time series considered in Mastroeni et al. (2018), the paper Rivero and Vellucci (2022) modeled the copper price time series with a matching pursuit algorithm, which used a waveform dictionary with rectangular window functions introduced in De Carli and Vellucci (2018). The novelty of the present paper is the introduction of an SDE for the copper futures market: to the best of our knowledge, this is the only paper in the literature that attempts to provide a mathematical model for the copper futures market.

This paper is organized as follows. Section 2 describes the state space reconstruction by introducing a formal definition of determinism for time series. In Section 3, we recall the whitening filter approach and prove some preliminary useful results. Section 4 contains a brief exploratory analysis of the data, whereas in Section 5, we give the main results of the paper, constructing an SDE characterized by a damped cosine autocorrelation function. In Section 6.2, we perform some numerical simulations on the SDE obtained in the preceding sections for copper futures prices’ time series and, in the Appendix A, we compare them with the results obtained from an artificial time series. Section 7 is devoted to our conclusions.

2. Background

2.1. State Space Reconstruction

With the state space reconstruction problem, we aim at recreating states where the only information available is contained in a time series. In other words, we would like the time series ${\left\{x\left(k\right)\right\}}_{k=1}^n=\{x\left(1\right),x\left(2\right),\dots ,x\left(n\right)\}$ to be approximately described by a smooth dynamical system f on a d-dimensional manifold $\mathcal{M}$ (which we assume to be $\mathcal{M}={\mathbb{R}}^d$):

$$s\left(t\right)=f^t\left(s\left(0\right)\right)$$

(1)

where $s\left(t\right)$ is the state at time t. For $t\ge 0$, $f^t$ is the time-t map of the continuous dynamical system, i.e., the transformation of state space which takes $s\left(0\right)$ to $s\left(t\right)$. For noise-free data, the time series is related to the dynamical system by

$$x\left(t\right)=h\left(s\right(t\left)\right),$$

(2)

where $h:\mathcal{M}\to \mathbb{R}$ is known as the measurement function.

To analyze a time series through the formalism of dynamical systems theory, the reconstruction of state space is necessarily the first step. There is no way to reconstruct states in their original form because f and h are typically both unknown but, luckily, the reconstructed state space is in some sense equivalent to the original.

The state space is reconstructed according to the delay coordinate method proposed by Takens (1981) and Packard et al. (1980). The delay coordinate method defines the state vector at time i as follows

$$\mathbf{x}\left(i\right)={\left[x\left(i\right),x(i+\tau ),x(i+2\tau ),\dots ,x(i+(d-1\left)\tau \right)\right]}^T,$$

(3)

where d is the embedding dimension, and $\tau$ is an appropriate time delay. Here, T denotes the transpose, and we adopt the convention that states are represented by column vectors. Obviously, if n is the length of the time series, then i ranges in $\{1,2,\cdots ,n-(d-1\left)\tau \right\}$. For a complete characterization of $\tau$ and d see Mastroeni et al. (2018, 2019); Mastroeni and Vellucci (2019, 2022) and references therein.

In Takens (1981), Takens studied the delay reconstruction map ${\mathsf{\Phi }}_{(f,h)}$, introducing the following result:

Theorem 1

(Takens 1981). For generic pairs $(f,h)$, where

• $f:\mathcal{M}\to \mathcal{M}$ is a $C^2$-diffeomorphism of $\mathcal{M}$ in itself, and
• $h:\mathcal{M}\to \mathbb{R}$ is a $C^2$-differentiable function,

the map ${\mathsf{\Phi }}_{(f,h)}:\mathcal{M}\to {\mathbb{R}}^{2d+1}$ defined by

$${\mathsf{\Phi }}_{(f,h)}\left(s\right)={\left[h\left(s\right),h\left(f^{\tau }\left(s\right)\right),\dots ,h\left(f^{2\tau d}\left(s\right)\right)\right]}^T$$

(4)

is an embedding of $\mathcal{M}$ in ${\mathbb{R}}^{2d+1}$.

The important point of the theorem depends on the topological notion of “genericity”:

Definition 1.

Let P be a property of functions in ${\mathcal{C}}^k(\mathcal{M},\mathcal{N})$, which they might or might not have. Then, we say that $P\left(f\right)$ is true for a generic $f\in C^k(\mathcal{M},\mathcal{N})$ if the set of functions for which it holds is open and dense in ${\mathcal{C}}^k$-topology.

For Theorem 1 to work, the genericity needs to be taken in the space of pairs $(f,h)$, but it can be replaced guaranteeing additional requirements for f, directed towards establishing a slightly different version of the theorem:

Theorem 2.

Let $f:\mathcal{M}\to \mathcal{M}$ be a $C^2$-diffeomorphism such that

• the periodic points with period $k\le 2d$ are finite in number;
• if x is any periodic point with period $k\le 2$ d, then the derivative of ${\varphi }^k$ at x has all distinct eigenvalues.

Then, for a generic $h\in C^2(\mathcal{M},\mathbb{R})$, the map ${\mathsf{\Phi }}_{(\varphi ,y)}:\mathcal{M}\to {\mathbb{R}}^{2d+1}$ defined as in Theorem 1 is an embedding of $\mathcal{M}$ in ${\mathbb{R}}^{2d+1}$.

Under the assumptions of Taken’s Theorem 1, an embedding is a smooth one-to-one coordinate transformation with a smooth inverse. For a compact manifold $\mathcal{M}$, a smooth one-to-one map $\mathcal{M}\to \mathcal{M}$ whose derivative is everywhere one-to-one is a smooth embedding, i.e., the inverse is smooth for free (the local smoothness of the inverse essentially comes from the inverse function theorem).

If $\mathsf{\Phi }$ is an embedding, then a smooth dynamics F is induced on the space of reconstructed vectors:

$$F^t\left(\mathbf{x}\right)=\mathsf{\Phi }\circ f^t\circ {\mathsf{\Phi }}^{-1}\left(\mathbf{x}\right),$$

(5)

where F is equivalent to the original dynamics f. Likewise, the reconstructed states can be used to estimate F; thus, we can use this for the purpose of the paper, as described in the following section.

2.2. A Definition of Determinism

In this section, we formalize the notion of determinism of a time series, taking it from the field of chaos theory. It is based on the recurrence plot, introduced by Eckmann et al. (1987) to visualize the different times at which a trajectory visits roughly the same area of the state space, which in turn is the result of the state space reconstruction theory, already seen in the previous section.

Let $\parallel \parallel$ be the Euclidean norm in ${\mathbb{R}}^d$ and $ϵ$ be a tolerance parameter to be chosen as recommended in Mastroeni et al. (2018, 2019); Mastroeni and Vellucci (2019, 2022) and references therein. Let

$$\mathbf{\Theta }\left(x\right):=\left\{\begin{array}{cc}1,\hfill & x>0\hfill \\ 0,\hfill & x\le 0\hfill \end{array}\right.$$

(6)

be the Heaviside step function. Then, we have the following definition of the recurrence plot.

Definition 2.

Let us consider the matrix

$$M_{ij}=\mathbf{\Theta }\left(ϵ-\parallel \mathit{v}\left(i\right)-\mathit{v}\left(j\right)\parallel \right).$$

(7)

The recurrence plot $\mathcal{R}$ is the set of all points $(i,j)$ for which $M_{ij}=1$:

$$\mathcal{R}:=\left\{(i,j)\in {\mathbb{R}}^2|M_{ij}=1\right\}$$

(8)

The matrix M is symmetric due to the symmetry of the Euclidean norm, so that the recurrence plot $\mathcal{R}$ is likewise symmetric around the bisecting line of the first quadrant.

Let $P\left(l\right)$ be the histogram of diagonal lines of length l, i.e.,

$$P\left(l\right)=\sum _{i,j=1}^N\left(1-M_{i-1,j-1}\right)\left(1-M_{i+l,j+l}\right)\prod _{k=0}^{l-1}{M}_{i+k,j+k}.$$

As recalled in Marwan et al. (2007), “processes with uncorrelated or weakly correlated, stochastic or chaotic behavior cause none or very short diagonals, whereas deterministic processes cause longer diagonals and less single, isolated recurrence points”. This allows us to introduce the following definition of the determinism of a system:

Definition 3.

Fix $l_{min}>1$. Then, the percentage of recurrence points placed on diagonal structures with respect to all recurrence points,

$$DET=\frac{{\sum }_{l=l_{min}}^{N}lP\left(l\right)}{{\sum }_{l=1}^{N}lP\left(l\right)},$$

(9)

is introduced as a measure for the determinism of the system.

(For the choice of $l_{min}$ see Marwan et al. (2007)). This is the measure for determinism of the reconstructed system and, from Equation (5), also of the smooth dynamical system f that approximately describes the starting time series.

3. Preliminary Results: The Whitening Filter Approach

Let us recall the whitening filter approach and introduce some preliminary results that we need to achieve the mathematical results in the subsequent sections.

One of the most commonly used methods to characterize a linear time-variant system, that is, to specify its input–output relationship, is based on the impulse response of the system. For example, a two-terminal pair system, with input $x\left(t\right)$ and output $y\left(t\right)$ is depicted in Figure 1. The impulse response of the linear time-variant system is denoted by $h(t,\tau )$. It is the output at time t to an input applied at time $\tau$. Given $T>0$, for a physically realizable system, we assume that $h(t,\tau )$ is zero outside the range $[0,T]$. Then, the input–output relation assumes the form:

$$y\left(t\right)={\int }_0^{T}h(t,\tau )x\left(\tau \right)d\tau$$

(10)

Let us consider now a nonstationary noise $n\left(t\right)$ that can be characterized by its autocorrelation $R_n(t,\tau )$. The whitening filter problem consists in finding the system $h(t,\tau )$ for a given $R_n(t,\tau )$, such that the output spectrum is $R_{n^{\prime }}(t-\tau )=\delta (t-\tau )$. Let us introduce the following Proposition.

Proposition 1.

Let $\left\{{\varphi }_j\right\}$ be the set of eigenfunctions of $R_n(\tau ,{\tau }^{\prime })$,

$${\lambda }_j{\varphi }_j\left(\tau \right)={\int }_0^T{R}_n(\tau ,{\tau }^{\prime }){\varphi }_j\left({\tau }^{\prime }\right)d{\tau }^{\prime }.$$

(11)

Then, the filter given by the expansion of equation:

$$h(t,\tau )=\sum _j{a}_j\left(t\right){\varphi }_j\left(\tau \right)where{a}_j\left(t\right)={\int }_0^{T}h(t,\tau ){\varphi }_j\left(\tau \right)d\tau$$

(12)

solves the whitening filter problem.

Proof. 

Given the expression of noise output

$$n^{\prime }\left(t\right)={\int }_0^{T}h(t,\tau )n\left(\tau \right)d\tau ,$$

(13)

we can compute its autocorrelation function as:

$$\begin{array}{cc}\hfill R_{n^{\prime }}(t_1,t_2)& =\mathbb{E}\left[n^{\prime }\left(t_1\right)n^{\prime }\left(t_2\right)\right]\hfill \\ \hfill & =\mathbb{E}\left[{\int }_0^{T}h(t_1,\tau )n\left(\tau \right)d\tau {\int }_0^{T}h(t_2,{\tau }^{\prime })n\left({\tau }^{\prime }\right)d{\tau }^{\prime }\right]\hfill \\ \hfill & =\mathbb{E}\left[{\int }_0^T{\int }_0^{T}h(t_1,\tau )h(t_2,{\tau }^{\prime })n\left(\tau \right)n\left({\tau }^{\prime }\right)d\tau d{\tau }^{\prime }\right]\hfill \\ \hfill & ={\int }_0^T{\int }_0^{T}h(t_1,\tau )h(t_2,{\tau }^{\prime })\mathbb{E}\left[n\left(\tau \right)n\left({\tau }^{\prime }\right)\right]d\tau d{\tau }^{\prime },\hfill \end{array}$$

(14)

and then,

$$R_{n^{\prime }}(t_1,t_2)={\int }_0^T{\int }_0^{T}h(t_1,\tau )h(t_2,{\tau }^{\prime })R_n(\tau ,{\tau }^{\prime })d\tau d{\tau }^{\prime }.$$

(15)

If $n\left(t\right)$ is real $R_n(\tau ,{\tau }^{\prime })$, it is symmetric; therefore, a standard eigenfunction approach can be used to solve the integral Equation (15).

Let $\left\{{\varphi }_j\right\}$ be the set of eigenfunctions of $R_n(\tau ,{\tau }^{\prime })$,

$${\lambda }_j{\varphi }_j\left(\tau \right)={\int }_0^T{R}_n(\tau ,{\tau }^{\prime }){\varphi }_j\left({\tau }^{\prime }\right)d{\tau }^{\prime }.$$

(16)

Since the set of eigenfunctions is complete, the system functions can be expanded as

$$h(t,\tau )=\sum _j{a}_j\left(t\right){\varphi }_j\left(\tau \right),$$

(17)

where

$$a_j\left(t\right)={\int }_0^{T}h(t,\tau ){\varphi }_j\left(\tau \right)d\tau .$$

(18)

By using the expansion (17), Equation (15) becomes:

$$\begin{array}{cc}\hfill R_{n^{\prime }}(t_1,t_2)& =\sum _j{a}_j\left(t_2\right){\int }_0^{T}h(t_1,\tau )\left[{\int }_0^T{R}_n(\tau ,{\tau }^{\prime }){\varphi }_j\left({\tau }^{\prime }\right)d{\tau }^{\prime }\right]d\tau \hfill \\ \hfill & =\sum _j{\lambda }_j{a}_j\left(t_2\right){\int }_0^{T}h(t_1,\tau ){\varphi }_j\left(\tau \right)d\tau \hfill \\ \hfill & =\sum _j{\lambda }_j{a}_j\left(t_1\right)a_j\left(t_2\right).\hfill \end{array}$$

(19)

If one first considers making the output noise spectrum only stationary:

$$R_{n^{\prime }}(t_1,t_2)=R(t_1-t_2),$$

(20)

then $a_j\left(t_1\right)$ is identified as ${\left({\gamma }_j/{\lambda }_j\right)}^{1/2}{\psi }_j\left(t_1\right)$, where ${\psi }_j\left(t_1\right)$ and ${\gamma }_j$ are the eigenfunction and eigenvalue associated with $R(t_1-t_2)$, respectively:

$${\gamma }_j{\psi }_j\left(t_1\right)={\int }_0^{T}R(t_1-t_2){\psi }_j\left(t_2\right)d{t}_2.$$

(21)

□

Now, we denote by $\mathcal{L}\left(f\left(s\right)\right)={\mathcal{L}}_f\left(s\right)$ or simply $F\left(s\right)$ the bilateral Laplace transform of a function $f:\mathbb{R}\mapsto \mathbb{R}$.

If a random (wide-sense stationary) process $n\left(t\right)$ passes through a time-invariant filter $h\left(\tau \right)$, such as the one depicted in Figure 1, then the autocorrelation function given by Equation (15) becomes

$$R_{n^{\prime }}\left(t\right)={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{R}_n(t-\tau +{\tau }^{\prime })h\left(\tau \right)h\left({\tau }^{\prime }\right)d\tau d{\tau }^{\prime }.$$

(22)

For white noise input, $R_n\left(t\right)=\delta \left(t\right)$, and

$$R_{n^{\prime }}\left(t\right)={\int }_{-\infty }^{\infty }h(\tau -t)h\left(\tau \right)d\tau .$$

(23)

The power spectrum of the output is the Laplace transform of $R_{n^{\prime }}\left(t\right)$:

$$S_{n^{\prime }}\left(s\right)={\int }_{-\infty }^{\infty }{R}_{n^{\prime }}\left(t\right)e^{-st}dt={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }h(\tau -t)h\left(\tau \right)e^{-st}dtd\tau .$$

(24)

By integrating first with respect to t:

$${\int }_{-\infty }^{\infty }h(\tau -t)e^{-st}dt=e^{-s\tau }{\int }_{-\infty }^{\infty }h\left(\gamma \right)e^{s\gamma }d\gamma =e^{-s\tau }H(-s),$$

(25)

and then with respect to $\tau$, we obtain:

$$S_{n^{\prime }}\left(s\right)=H(-s){\int }_{-\infty }^{\infty }h\left(\tau \right)e^{-s\tau }d\tau =H(-s)H\left(s\right).$$

(26)

Let us now explain what the relation is between a pole or zero of $H\left(s\right)$ and the corresponding pole or zero of $H\left(-s\right)$. The poles and zeros of $H\left(-s\right)$ are those of $H\left(s\right)$ mirrored across the imaginary axis of a complex plane. Therefore, the poles and zeros of $S_{n^{\prime }}\left(s\right)$ are paired across the imaginary axis. With N pole pairs and M zero pairs in $S_{n^{\prime }}\left(s\right)$, we have (in principle) $2^{N+M}$ different choices of whitening filter that produce $S_{n^{\prime }}\left(s\right)$. However, imposing the causality/stability of the whitening filter, we have to place all the poles in the left-hand-side complex plane. After that, we have just $2^M$ choices. Then, to impose a minimum phase, we place all zeros in the left-hand-side complex plane. At the end, the filter is stable, causal, and minimum phase. Here, we refer to a strictly stable linear system as being one with all the poles of its transfer function in the open left-hand-side complex plane ${\mathbb{C}}_{0^{-}}=\{s\in \mathbb{C}:\Re s<0\}$, and we refer to a stable linear system as being one with poles in the closed left-hand-side complex plane.

4. Exploratory Data Analysis

The dataset employed in our analysis was the time series of intraday Generic 1st Futures Copper closing prices (HG1 ticker) as exchanged on the COMEX market (CMX) and retrieved from the Bloomberg website. The dataset consisted of 2864 observations and spanned from 8 June 2021 to 8 June 2022.

We denoted the time series of intraday prices by ${\left\{P\left(t\right)\right\}}_{t=1}^n$, with $P\left(t\right)\in {\mathbb{R}}^{+}$, and the logarithmic returns by ${\left\{p\left(t\right)\right\}}_{t=1}^{n-1}={\left\{ln\frac{P(t+1)}{P\left(t\right)}\right\}}_{t=1}^{n-1}$. The overall time series is shown in Figure 2 and Figure 3.

The data were not equally spaced intraday prices, i.e., the number of intraday observations was not constant but varied over time. See Figure 4.

The Hurst exponent of the series $p\left(t\right)$ was 0.53, which meant that the time series was trending. This value was near to $0.5$ (i.e., the time series was generated by a geometric Brownian motion).

Figure 5 shows the autocorrelation of the considered data set for time lags up to points in time 10 (number of observations, recall that we have intraday prices). The red line is the autocorrelation computed from the data, while the green line corresponds to the nonlinear fit according to the expression

$$R\left(\tau \right)={\sigma }^2{e}^{-\alpha \left|\tau \right|}cos\left(\omega \tau \right)$$

(27)

The damped cosine function was a good approximation to the autocorrelation for time lags up to approximately 10 observations.

5. Construction of the SDE-Based Model

Let us assume the logarithmic returns $p\left(t\right)=ln\frac{P(t+1)}{P\left(t\right)}$ were characterized by the damped cosine autocorrelation $R\left(\tau \right)=\mathbb{E}\left[p\left(t\right)p(t+\tau )\right]$ described by (27), whose decay rate is governed by the coefficient $\alpha$.

Since we know that the autocorrelation price process $p\left(t\right)$ fits a function such as (27), the problem was to derive an SDE that describes a process with this autocorrelation. Theorem 3 answers this question:

Theorem 3.

Let $w\left(t\right)$ be the white noise. Let also $p\left(0\right)={\kappa }_1$, $\dot{p}\left(0\right)={\kappa }_2$, and $w\left(0\right)={\kappa }_3$. If $R_p\left(\tau \right)=\mathbb{E}\left[p\left(t\right)p(t+\tau )\right]$ is the autocorrelation of log-price $p\left(t\right)$ described by (27), then the SDE that results from the whitening filter procedure is

$$\ddot{p}\left(t\right)+2\alpha \dot{p}\left(t\right)+{\beta }^{2}p\left(t\right)+\left({\kappa }_2-{\kappa }_3+2\alpha {\kappa }_1\right)\delta \left(t\right)+{\kappa }_1\dot{\delta }\left(t\right)=\beta w\left(t\right)+\dot{w}\left(t\right),$$

(28)

where δ is the Dirac delta function, $\beta =\sqrt{{\omega }^2+{\alpha }^2}$, and the variance of $w\left(t\right)$ is $2\alpha {\sigma }^2$.

Proof. 

The Wiener–Kolmogorov whitening procedure from Wiener (1949) is used on $R_p\left(\tau \right)$ to represent the log-price $p\left(t\right)=ln\frac{P(t+1)}{P\left(t\right)}$ as a function of time, driven by a white noise input $w\left(t\right)$. The procedure is used in the proof to decompose the Laplace transform of $R_p\left(\tau \right)$ into the product of white noise and a system function: ${\mathcal{L}}_R\left(s\right)=W\left(s\right)H(-s)H\left(s\right)$. Here, the quantity $H\left(s\right)$ is the transfer function of the whitening filter for the log-price $p\left(t\right)$, and $W\left(s\right)$ is the transform of the white noise $w\left(t\right)$ that drives $p\left(t\right)$.

The Laplace transform of the autocorrelation function $R_p\left(\tau \right)$ is

$${\mathcal{L}}_R\left(s\right)={\int }_{-\infty }^{\infty }{\sigma }^2{e}^{-\alpha \left|\tau \right|}cos\left(\omega \tau \right)e^{-s\tau }d\tau .$$

(29)

Splitting the module,

$${\mathcal{L}}_R\left(s\right)={\int }_0^{\infty }{\sigma }^{2}cos\left(\omega \tau \right)e^{-(\alpha -s)\tau }d\tau +{\int }_0^{\infty }{\sigma }^{2}cos\left(\omega \tau \right)e^{-(\alpha +s)\tau }d\tau ,$$

(30)

solving the integrals,

$${\mathcal{L}}_R\left(s\right)={\sigma }^2\left(\frac{\alpha -s}{{\omega }^2+{(\alpha -s)}^2}+\frac{s+\alpha }{{\omega }^2+{(s+\alpha )}^2}\right),$$

(31)

and factorizing the resulting expression, we obtain

$${\mathcal{L}}_R\left(s\right)=2\alpha {\sigma }^2\left(\frac{\beta -s}{s^2-2\alpha s+{\beta }^2}\right)\left(\frac{\beta +s}{s^2+2\alpha s+{\beta }^2}\right),$$

(32)

where $\beta =\sqrt{{\omega }^2+{\alpha }^2}$. Hence, we have

$$H\left(s\right)=\frac{\beta +s}{s^2+2\alpha s+{\beta }^2},$$

(33)

and

$$W\left(s\right)=2\alpha {\sigma }^2.$$

(34)

Indeed, the autocorrelation function $R_w\left(\tau \right)$ of the white noise input is found by taking the inverse Laplace transform of (34) as

$$R_w\left(\tau \right)=2\alpha {\sigma }^2\delta \left(\tau \right),$$

(35)

where $2\alpha {\sigma }^2$ is the variance of $w\left(t\right)$.

From the definition of a transfer function, we have that $H\left(s\right)=\frac{P\left(s\right)}{W\left(s\right)}$, where $P\left(s\right)$ is the Laplace transform of $p\left(t\right)$. Moreover, it is well known that for a generic function $y\left(t\right)$,

$$\begin{array}{cc}\hfill sY\left(s\right)& =\mathcal{L}\left(\dot{y}\left(t\right)\right)+y\left(0\right)=\mathcal{L}(\dot{y}\left(t\right)+y\left(0\right)\delta \left(t\right))\hfill \\ \hfill s^{2}Y\left(s\right)& =\mathcal{L}\left(\ddot{y}\left(t\right)\right)+sy\left(0\right)+\dot{y}\left(0\right)=\mathcal{L}(\ddot{y}\left(t\right)+y\left(0\right)\dot{\delta }\left(t\right)+\dot{y}\left(0\right)\delta \left(t\right)).\hfill \end{array}$$

Applying these last considerations to Equation (33), and recalling that $p\left(0\right)={\kappa }_1$, $\dot{p}\left(0\right)={\kappa }_2$, and $w\left(0\right)={\kappa }_3$, we obtain the differential Equation (28). □

Remark 1.

As anticipated in Section 3, it can be easily demonstrated that the whitening filter, whose transfer function is described by Equation (33) in the proof of Theorem 3, is unique, stable, causal, and minimum phase.

At this stage, we do not yet have the formulation of a standard SDE, i.e., a rigorously defined equivalent implicit Ito’s integral equation. For instance,

$$\begin{array}{cc}\hfill dp\left(t\right)& =a\left(t,p\left(t\right)\right)dt+b\left(t,p\left(t\right)\right)dB\left(t\right)\iff \hfill \\ \hfill p\left(t\right)& =p\left(0\right)+{\int }_0^{t}a\left(s,p\left(s\right)\right)ds+{\int }_0^{t}b\left(s,p\left(s\right)\right)d{B}_s,\hfill \end{array}$$

(36)

where B denotes a Wiener process (the standard Brownian motion). Instead, the approach followed in Theorem 3 presents us with these mathematical problems Evans (2012):

• Defining the white noise $w\left(t\right)$ in a rigorous way;
• Showing that (28) has a solution, discussing uniqueness, asymptotic behavior, and dependence upon $\alpha$ and $\omega$.

However, we can start from (28) with a heuristic explanation: we could think of the white noise introduced in Theorem 3 as the “derivative” of Brownian motion. It is just heuristic reasoning because the standard Brownian motion is nowhere differentiable with probability one, but it is a way to start and move intuitively from (28). To formalize this rigorously, we use the theory of generalized functions.

As we know, a generalized function is a continuous linear function on the set of infinitely differentiable functions with bounded support, which we denote with $C_0^{\infty }$ or simply $\mathcal{D}$. The space of generalized functions on $\mathcal{D}$ is denoted by ${\mathcal{D}}^{\prime }$. We write

$$F\left[\varphi \right]:\mathcal{D}\to \mathbb{R}.$$

(37)

We define the generalized derivative of $F\left[\varphi \right]$ to be the generalized function $F^{\prime }\left[\varphi \right]$, so that

$$F^{\prime }\left[\varphi \right]\equiv -F\left[{\varphi }^{\prime }\right],\mathrm{for}\mathrm{all}\varphi \in \mathcal{D}.$$

(38)

Thinking of Brownian motion as a random function, we can define a white noise $w\left(t\right)$ as its generalized derivative. Hence:

Lemma 1.

Let $B\left(t\right)\ge 0$ be a real-valued Brownian motion on a probability space $(\mathbf{\Omega },\mathcal{A},P)$. Then,

$$\dot{B}\left(t\right)=w\left(t\right).$$

(39)

Proof. 

Formally, according to Evans (2012); Lord et al. (2014), white noise $w\left(t\right)$ is a stochastic process on a probability space $(\mathbf{\Omega },\mathcal{A},P)$ with $\mathbb{E}\left[w\right(t\left)\right]=0$ and $\mathbb{E}\left[w\right(s\left)w\right(t\left)\right]=\delta (s-t)$ for all $s,t\in [0,1]$, where $\delta$ denotes the Dirac delta function.

We define

$$\mathrm{D}(t,\Delta t):=\frac{B(t+\Delta t)-B\left(t\right)}{\Delta t}\mathrm{for}t\ge 0,$$

(40)

for some $\Delta t>0$. Then, $\mathrm{D}(t,\Delta t)$ is a stochastic process on $(\mathbf{\Omega },\mathcal{A},P)$ with

$$\mathbb{E}\left[\mathrm{D}(t,\Delta t)\right]=0\mathrm{for}\mathrm{all}t\ge 0,$$

(41)

and for all $s,t\ge 0$,

$$\eta (s-t,\Delta t):=Cov\left[\mathrm{D}(s,\Delta t),\mathrm{D}(t,\Delta t)\right]=\left\{\begin{array}{cc}{\displaystyle \frac{\Delta t-|s-t|}{\Delta t^2}}\hfill & ,\mathrm{if}|s-t|\le \Delta t\hfill \\ 0\hfill & ,\mathrm{if}|s-t|\ge \Delta t\hfill \end{array}\right..$$

(42)

Let $\varphi \in \mathcal{D}$. Then, we have

$$\underset{\Delta t\to 0}{lim}{\int }_{-\infty }^{\infty }\varphi \left(x\right)\eta (x,\Delta t)dx=\underset{\Delta t\to 0}{lim}{\int }_{-\Delta t}^{\Delta t}\varphi \left(x\right)\frac{\Delta t-\left|x\right|}{\Delta t^2}dx,$$

(43)

i.e., replacing $x\mapsto \Delta tx$,

$$\underset{\Delta t\to 0}{lim}{\int }_{-1}^1\varphi (\Delta tx)(1-|x\left|\right)dx=\varphi \left(0\right){\int }_{-1}^1(1-|x\left|\right)dx=\varphi \left(0\right).$$

(44)

Therefore, in distribution, we can write

$$\eta (x,\Delta t)\stackrel{\Delta t\to 0}{\to }\delta \left(x\right).$$

(45)

Thus, we have that the derivative of the covariance matrix for the Brownian motion is the same as the covariance matrix for the white noise; thus, the derivative of the Brownian motion is the white noise. □

In a similar manner, we introduce the second time derivative on a Wiener process (or just the time derivative on a white noise):

Lemma 2.

Let $B\left(t\right)\ge 0$ be a real-valued Brownian motion on a probability space $(\mathbf{\Omega },\mathcal{A},P)$. Then, the stochastic process

$$\ddot{B}\left(t\right)=\dot{w}\left(t\right)=\eta \left(t\right)$$

(46)

has the following statistics:

$$\mathbb{E}\left[\eta \left(t\right)\right]=0;\mathbb{E}\left[\eta \left(t\right)\eta \left(s\right)\right]=-\frac{d^2}{d{z}^2}\delta \left(z\right){|}_{z=t-s}.$$

(47)

Proof. 

Similar to the proof of Lemma 1, we define

$$\mathrm{D}(t,\Delta t):=\frac{w(t+\Delta t)-w\left(t\right)}{\Delta t}\mathrm{for}t\ge 0,$$

(48)

and we take $\Delta t$ at the last step.

For the expected value of $\mathrm{D}(t,\Delta t)$ we have that $\frac{1}{\Delta t}\left\{\mathbb{E}\left[w(t+\Delta t)\right]-\mathbb{E}\left[w\left(t\right)\right]\right\}=0$. For the covariance, we have instead

$$\mathbb{E}\left[\mathrm{D}(t,\Delta t)\mathrm{D}(s,\Delta t)\right]=\frac{1}{\Delta t^2}\mathbb{E}\left[\left(w(t+\Delta t)-w\left(t\right)\right)\left(w(s+\Delta t)-w\left(s\right)\right)\right],$$

(49)

i.e.,

$$\begin{array}{cc}\hfill \mathbb{E}\left[\mathrm{D}(t,\Delta t)\mathrm{D}(s,\Delta t)\right]& =\frac{1}{\Delta t^2}\{\mathbb{E}\left[w(t+\Delta t)w(s+\Delta t)\right]-\mathbb{E}\left[w(t+\Delta t)w\left(s\right)\right]+\hfill \\ \hfill & -\mathbb{E}\left[w\right(t\left)w\right(s+\Delta t\left)\right]+\mathbb{E}\left[w\right(t\left)w\right(s\left)\right]\}.\hfill \end{array}$$

(50)

From the statistics of a white noise,

$$\begin{array}{cc}\hfill \mathbb{E}\left[w\right(t\left)w\right(s\left)\right]& =\mathbb{E}\left[w\right(t+\Delta t\left)w\right(s+\Delta t\left)\right]=\delta (t-s)\hfill \\ \hfill \mathbb{E}\left[w\right(t+\Delta t\left)w\right(s\left)\right]& =\delta (t-s+\Delta t)\hfill \\ \hfill \mathbb{E}\left[w\right(t\left)w\right(s+\Delta t\left)\right]& =\delta (t-s-\Delta t),\hfill \end{array}$$

(51)

we have that

$$\begin{array}{cc}\hfill \mathbb{E}\left[\mathrm{D}(t,\Delta t)\mathrm{D}(s,\Delta t)\right]& =-\frac{1}{\Delta t^2}\left[\delta (t-s+\Delta t)-2\delta (t-s)+\delta (t-s-\Delta t)\right]\hfill \end{array}$$

(52)

$$\begin{array}{cc}\hfill & =-\frac{d^2}{d{z}^2}\delta \left(z\right){|}_{z=t-s}\hfill \end{array}$$

(53)

as $\Delta t\to 0$. □

Corollary 1.

The Ito’s integral form of Equation (28) is

$$\left\{\begin{array}{c}\mathrm{d}p\left(t\right)=v\left(t\right)\mathrm{d}t-{\kappa }_1\mathrm{d}\mathbf{\Theta }+\mathrm{d}B\left(t\right),\left(p\left(0\right)={\kappa }_1,v\left(0\right)=2{\kappa }_1-{\kappa }_3\right)\hfill \\ \mathrm{d}v\left(t\right)=-\left({\beta }^{2}p\left(t\right)+2\alpha v\left(t\right)\right)\mathrm{d}t-\left({\kappa }_2-{\kappa }_3\right)\mathrm{d}\mathbf{\Theta }+(\beta -2\alpha )\mathrm{d}B\left(t\right)\hfill \end{array}\right.,$$

(54)

where Θ is the Heaviside step function with the distributional derivative $\dot{\mathbf{\Theta }}=\delta$. Equation (54) has the unique solution:

$$\left\{\begin{array}{c}p\left(t\right)=\frac{e^{-\alpha t}}{\omega }\left[p_0\left(\omega cos\omega t+\alpha sin\omega t\right)+v_{0}sin\omega t\right]+\hfill \\ +\frac{e^{-\alpha t}}{\omega }{\int }_0^t{e}^{\alpha s}\left(\omega cos\omega (t-s)+(\beta -\alpha )sin\omega (t-s)\right)\mathrm{d}B\left(s\right)+\hfill \\ +\frac{e^{-\alpha t}}{\omega }\left(-{\kappa }_1\left(\omega cos\omega t+\alpha sin\omega t\right)+({\kappa }_3-{\kappa }_2)sin\omega t\right)\hfill \\ v\left(t\right)=\frac{e^{-\alpha t}}{\omega }\left[p_0{\beta }^{2}sin\omega t+v_0\left(\omega cos\omega t-\alpha sin\omega t\right)\right]+\hfill \\ +\frac{e^{-\alpha t}}{\omega }{\int }_0^t{e}^{\alpha s}\left(\omega (\beta -2\alpha )cos\omega (t-s)+({\beta }^2-\alpha \beta +2{\alpha }^2)sin\omega (t-s)\right)\mathrm{d}B\left(s\right)\hfill \\ +\frac{e^{-\alpha t}}{\omega }\left(-{\kappa }_1{\beta }^{2}sin\omega t+({\kappa }_3-{\kappa }_2)\left(\omega cos\omega t-\alpha sin\omega t\right)\right)\hfill \end{array}\right..$$

(55)

Proof. 

From Lemmas 1 and 2, we can rewrite Equation (28) as

$$\ddot{p}\left(t\right)+2\alpha \dot{p}\left(t\right)+{\beta }^{2}p\left(t\right)+\left({\kappa }_2-{\kappa }_3+2\alpha {\kappa }_1\right)\delta \left(t\right)+{\kappa }_1\dot{\delta }\left(t\right)=\beta \dot{B}\left(t\right)+\ddot{B}\left(t\right).$$

(56)

To reduce the order in Equation (56), we combine the derivative terms shifting some of them around. Hence, we obtain a first-order system with the substitution $v\left(t\right)=\dot{p}\left(t\right)-\dot{B}\left(t\right)+{\kappa }_1\delta \left(t\right)$, i.e., $\dot{v}\left(t\right)+\ddot{B}\left(t\right)-{\kappa }_1\dot{\delta }\left(t\right)=\ddot{p}$. Accordingly,

$$\begin{array}{cc}\hfill \dot{v}\left(t\right)+\ddot{B}\left(t\right)-{\kappa }_1\dot{\delta }\left(t\right)& +2\alpha \left(v\left(t\right)+\dot{B}\left(t\right)-{\kappa }_1\delta \left(t\right)\right)+{\beta }^{2}p\left(t\right)+\hfill \\ \hfill & +\left({\kappa }_2-{\kappa }_3+2\alpha {\kappa }_1\right)\delta \left(t\right)+{\kappa }_1\dot{\delta }\left(t\right)=\beta \dot{B}\left(t\right)+\ddot{B}\left(t\right)\Rightarrow .\hfill \end{array}$$

(57)

With simple algebra, we obtain

$$\left\{\begin{array}{c}\dot{p}\left(t\right)=v\left(t\right)-{\kappa }_1\delta \left(t\right)+\dot{B}\left(t\right)\hfill \\ \dot{v}\left(t\right)=-2\alpha v\left(t\right)-{\beta }^{2}p\left(t\right)-\left({\kappa }_2-{\kappa }_3\right)\delta \left(t\right)+\left(\beta -2\alpha \right)\dot{B}\left(t\right)\hfill \end{array}\right.,$$

(58)

and thus the corresponding well-defined Ito’s integral form (54).

By passing to the vectorial form and neglecting the Heaviside terms, Equation (54) can be rewritten as

$$d\left[\begin{array}{c}p\left(t\right)\hfill \\ v\left(t\right)\hfill \end{array}\right]=A\left[\begin{array}{c}p\left(t\right)\hfill \\ v\left(t\right)\hfill \end{array}\right]dt+\left[\begin{array}{c}1\hfill \\ \beta -2\alpha \hfill \end{array}\right]dB\left(t\right),$$

(59)

where

$$A=\left[\begin{array}{cc}0& 1\\ -{\beta }^2& -2\alpha \end{array}\right].$$

(60)

Given any initial value $(p\left(0\right),v\left(0\right))=\left(p_0,v_0\right)\in {\mathbb{R}}^2$, from the linearity of the SDE, we know that Equation (59) has the unique solution Mao (2007)

$$\left[\begin{array}{c}p\left(t\right)\hfill \\ v\left(t\right)\hfill \end{array}\right]=e^{At}\left[\begin{array}{c}{p}_0\hfill \\ v_0\hfill \end{array}\right]+{\int }_0^t{e}^{A(t-s)}\left[\begin{array}{c}1\hfill \\ \beta -2\alpha \hfill \end{array}\right]dB\left(s\right).$$

(61)

To consider the neglected terms, let us replace the Brownian motion $B\left(t\right)$ by the process $X\left(t\right)$ such that

$$\begin{array}{cc}\hfill X\left(t\right):& ={\int }_0^t\left[\begin{array}{c}1\hfill \\ \beta -2\alpha \hfill \end{array}\right]dB\left(s\right)+{\int }_0^t\left[\begin{array}{c}-{\kappa }_1\hfill \\ {\kappa }_3-{\kappa }_2\hfill \end{array}\right]d\mathbf{\Theta }\left(s\right)\hfill \\ \hfill & =\left[\begin{array}{c}1\hfill \\ \beta -2\alpha \hfill \end{array}\right]B\left(t\right)+\left[\begin{array}{c}-{\kappa }_1\hfill \\ {\kappa }_3-{\kappa }_2\hfill \end{array}\right]\mathbf{\Theta }\left(t\right)+X\left(0\right)=\left[\begin{array}{c}{X}^{\left(1\right)}\left(t\right)\hfill \\ X^{\left(2\right)}\left(t\right)\hfill \end{array}\right],\hfill \end{array}$$

(62)

whose quadratic variation is

$$\begin{array}{cc}\hfill {\left[X^{\left(1\right)},X^{\left(1\right)}\right]}_t& =t+{\kappa }_1^2\mathbf{\Theta }\left(t\right),{\left[X^{\left(2\right)},X^{\left(2\right)}\right]}_t={(\beta -2\alpha )}^{2}t+{({\kappa }_3-{\kappa }_2)}^2\mathbf{\Theta }\left(t\right)\hfill \\ \hfill {\left[X^{\left(1\right)},X^{\left(2\right)}\right]}_t& =(\beta -2\alpha )t+{\kappa }_1({\kappa }_2-{\kappa }_3)\mathbf{\Theta }\left(t\right).\hfill \end{array}$$

(63)

Now, for a one-dimensional Brownian motion $\left\{B\right(t):t\ge 0\}$, we have

$$\begin{array}{cc}\hfill \mathbb{E}\left[B\left(t\right)\mid {\mathcal{F}}^{+}\left(s\right)\right]& =\mathbb{E}\left[B\left(t\right)-B\left(s\right)\mid {\mathcal{F}}^{+}\left(s\right)\right]+B\left(s\right)\hfill \end{array}$$

(64)

$$\begin{array}{cc}\hfill & =\mathbb{E}\left[B\right(t)-B(s\left)\right]+B\left(s\right)=B\left(s\right).\hfill \end{array}$$

(65)

Hence, a Brownian motion is a martingale with respect to the filtration $\left({\mathcal{F}}^{+}\left(t\right):t\ge 0\right)$. Moreover, using the definition of the total variation of a real-valued function, it is possible to show that the variation of $\mathbf{\Theta }$ is 0 on every interval in $(-\infty ,0)$ and 1 on every larger interval. This means that the total variation of $\mathbf{\Theta }$ is $\mathbf{\Theta }$ itself; hence, it is a càdlàg adapted process of finite variation. Thus, $X\left(t\right)$ is a semi martingale Rogers and Williams (2000).

Hence, Equation (54) can be rewritten as

$$d\left[\begin{array}{c}p\left(t\right)\hfill \\ v\left(t\right)\hfill \end{array}\right]=A\left[\begin{array}{c}p\left(t\right)\hfill \\ v\left(t\right)\hfill \end{array}\right]dt+dX\left(t\right).$$

(66)

Ito’s formula gives us the solution of (66), i.e.,

$$\left[\begin{array}{c}p\left(t\right)\\ v\left(t\right)\end{array}\right]=e^{At}\left[\begin{array}{c}{p}_0\\ v_0\end{array}\right]+{\int }_0^t{e}^{A(t-s)}dX\left(s\right).$$

(67)

The matrix A has the property that $A\left[\begin{array}{c}p\\ v\end{array}\right]=\left[\begin{array}{c}v\\ -{\beta }^{2}p-2\alpha v\end{array}\right],$ so we can use it to define sequences $x_n,y_n$ as

$$\begin{array}{cc}\hfill & A^n\left[\begin{array}{c}1\\ 0\end{array}\right]=\left[\begin{array}{c}{x}_n\\ x_{n+1}\end{array}\right]\hfill \\ \hfill & A^n\left[\begin{array}{c}0\\ 1\end{array}\right]=\left[\begin{array}{c}{y}_n\\ y_{n+1}\end{array}\right],\hfill \end{array}$$

where $x_0=1,x_1=0,y_0=0,y_1=1$, and

$$\begin{array}{cc}\hfill & x_{n+2}=-2\alpha x_{n+1}-{\beta }^2{x}_n\hfill \\ \hfill & y_{n+2}=-2\alpha y_{n+1}-{\beta }^2{y}_n.\hfill \end{array}$$

The general closed formula for $x_n,y_n$, when the discriminant of the quadratic polynomial $\Delta$ is not zero, is in terms of $x_{+},x_{-}=\frac{-2\alpha \pm \sqrt{\Delta }}{2}$. In our case, $\Delta =4{\alpha }^2-4{\beta }^2=-4{\omega }^2$, which is by definition different from 0. Thus, $x_{+},x_{-}=-\alpha \pm i\omega$, and

$$\begin{array}{cc}\hfill & y_n=\frac{1}{2\omega i}\left(x_{+}^n-x_{-}^n\right)\hfill \\ \hfill & x_n=\frac{1}{2\omega i}\left(x_{+}{x}_{-}^n-x_{-}{x}_{+}^n\right).\hfill \end{array}$$

So:

$$\sum _{k=0}^{\infty }\frac{1}{k!}{\left(At\right)}^k\left[\begin{array}{c}0\\ 1\end{array}\right]=\left[\begin{array}{c}\sum _{k=0}^{\infty }\frac{y_k{t}^k}{k!}\\ \sum _{k=0}^{\infty }\frac{y_{k+1}{t}^k}{k!}\end{array}\right]=\frac{1}{2\omega i}\left[\begin{array}{c}{e}^{x_{+}t}-e^{x_{-}t}\\ x_{+}{e}^{x_{+}t}-x_{-}{e}^{x_{-}t}\end{array}\right],$$

which is the right column of $e^{At}.$ The left column is similarly calculated as:

$$\frac{1}{2\omega i}\left[\begin{array}{c}{x}_{+}{e}^{x_{-}t}-x_{-}{e}^{x_{+}t}\\ -{\beta }^2\left(e^{x_{-}t}-e^{x_{+}t}\right)\end{array}\right].$$

Equation (54) then has the unique solution:

$$\left\{\begin{array}{c}p\left(t\right)=\frac{e^{-\alpha t}}{\omega }\left[p_0\left(\omega cos\omega t+\alpha sin\omega t\right)+v_{0}sin\omega t\right]+\hfill \\ +\frac{e^{-\alpha t}}{\omega }{\int }_0^t{e}^{\alpha s}\left(\omega cos\omega (t-s)+(\beta -\alpha )sin\omega (t-s)\right)\mathrm{d}B\left(s\right)+\hfill \\ +\frac{e^{-\alpha t}}{\omega }{\int }_0^t{e}^{\alpha s}\left(-{\kappa }_1\left(\omega cos\omega (t-s)+\alpha sin\omega (t-s)\right)+({\kappa }_3-{\kappa }_2)sin\omega (t-s)\right)\mathrm{d}\mathbf{\Theta }\left(s\right)\hfill \\ v\left(t\right)=\frac{e^{-\alpha t}}{\omega }\left[p_0{\beta }^{2}sin\omega t+v_0\left(\omega cos\omega t-\alpha sin\omega t\right)\right]+\hfill \\ +\frac{e^{-\alpha t}}{\omega }{\int }_0^t{e}^{\alpha s}\left(\omega (\beta -2\alpha )cos\omega (t-s)+({\beta }^2-\alpha \beta +2{\alpha }^2)sin\omega (t-s)\right)\mathrm{d}B\left(s\right)\hfill \\ +\frac{e^{-\alpha t}}{\omega }{\int }_0^t{e}^{\alpha s}\left(-{\kappa }_1{\beta }^{2}sin\omega (t-s)+({\kappa }_3-{\kappa }_2)\left(\omega cos\omega (t-s)-\alpha sin\omega (t-s)\right)\right)\mathrm{d}\mathbf{\Theta }\left(s\right)\hfill \end{array}\right..$$

(68)

Thus, since $\mathbf{\Theta }\left(0\right)=0$, we have the claim. □

6. Numerical Implementation and Experiments

In this section, we demonstrate how the whitening filter approach with a damped cosine function as an autocorrelation function was able to fit the data through examples of simple simulations and intraday prices of copper futures.

6.1. Simulations

We generated an artificial time series that was a damped cosine function, whose oscillations were blurred by white noise $w\left(t\right)$. The time series was the following:

$$s\left(t\right)=2e^{-\frac{t}{100}}cos\left(\frac{t}{2}\right)w\left(t\right),$$

(69)

where $w\left(t\right)\in [0,1]$$\forall t$. Figure 6 shows the autocorrelation of $s\left(t\right)$ for time lags up to points in time 30. The red line is the autocorrelation computed from data, while the green line corresponds to a nonlinear fit, according to expression (27). Compared to Figure 5 for copper prices, the damped cosine function perfectly fit the empirical autocorrelation of the artificial data.

As shown in Figure 7, we randomly generated 250 observations from the distribution (69), and our approximation (green), according to model (54), was able to come close to the true values of the series up to about 60 observations. This was certainly a better result than that observed for the copper price series. The improved results obtained on the artificial series were due to the fact that it had an empirical autocorrelation function that perfectly fit the damped cosine assumed by our model. Moreover, let us recall that, from Section 6.2, the determinism coefficient seemed to be able to explain the goodness of the fit. Indeed, the determinism of the artificial data was much higher than that of copper prices—DET $=0.88$ for $s\left(t\right)$.

6.2. Intraday Copper Futures Prices

First, time-series data were segmented into consecutive equal segments by applying a time window of width $n=250$ observations. All the segments were disjoint. In particular, we focused on the following time ranges: (i) from 8 June 2021 to 8 July 2021; (ii) from 9 July 2021 to 10 August 2021; (iii) from 21 September 2021 to 17 November 2021; (iv) from 18 November 2021 to 8 December 2021; (v) from 9 December 2021 to 12 January 2022; (vi) from 13 January 2022 to 28 February 2022; and (vii) from 28 April 2022 to 18 May 2022. All these ranges were characterized by an autocorrelation function, such as the one depicted in Figure 5. In the following, we denote them by the corresponding Roman numerals.

In order to perform the simulations for the SDE model (54), we employed the R solver deSolve. This package uses a scheme introduced by Petzold (1983) for automatically determining whether an initial value problem, $dy/dt=f(y,t),y\left(t_0\right)=y_0$, can be solved more efficiently using a class of methods suited for nonstiff problems or a class of methods designed for stiff problems. We approximated a Dirac delta by a Gaussian function whose width was controlled by a parameter $\sigma$. We chose a Gaussian function located in the origin and with $\sigma =0.1$. We also used a geometric Brownian motion simulation model with 2000 possible trajectories. We fixed one of these trajectories, adding it to the deterministic part with Dirac delta term, according to the scheme (54). We repeated this procedure for all the trajectories, selecting the one which minimized the sum of the residuals

$$r\left(\widehat{p}\left(t\right)\right):=\sum _{t=1}^{n-1}{\left(p\left(t\right)-\widehat{p}\left(t\right)\right)}^2,$$

(70)

where n is the width of the time window, $\left\{p\right(t\left)\right\}$ is the time series of logarithmic returns, and $\left\{\widehat{p}\left(t\right)\right\}$ represents the simulated prices according to Equation (54). We denote by ${\widehat{p}}_{\mathrm{opt}}\left(t\right)$ the simulated price series, which minimizes (70), and $r_{\mathrm{opt}}=r\left({\widehat{p}}_{\mathrm{opt}}\left(t\right)\right)$.

However, since the log returns of copper were cosine-exponentially autocorrelated only for time lags up to 10 observations, the ability of the developed model (54) to represent the copper price behavior was expected to be acceptable only for that time frame. Indeed this was the case, as shown in Figure 8, Figure 9 and Figure 10. The trajectories were very close to the real data up to 10 observations. The residuals estimated for the optimum trajectory, $r_{\mathrm{opt}}$, can be read in

Table 1. Characterization of the log returns of copper for different time ranges.

Range	Standard Deviation	Determinism	${\mathit{r}}_{\mathbf{opt}}$
(i)	0.00002467958	0.2163743	14.96101
(ii)	0.00001708622	0	18.56718
(iii)	0.00004668479	0.3684211	18.03019
(iv)	0.00001576487	0.4649123	15.72518
(v)	0.00001830713	0	21.52263
(vi)	0.00002218824	0.4850746	17.27463
(vii)	0.00002092907	0	28.45987

where we also detailed the values of the standard deviation and determinism for each of the considered time ranges. The model’s ability to fit data only up to 10 observations was related to the following points: (1) the data did not really show an inflection and sign change in the first three points, as we would expect from the cosine—see Figure 5—and the damped cosine function (27) was just a valid approximation for the first few lags; (2) the determinism level in Table 1 revealed the presence of a large stochastic component, and this meant that the behavior of the copper price could be predicted only in the short term Mastroeni and Vellucci (2019, 2022). Regarding Table 1, the worst values of $r_{\mathrm{opt}}$ were achieved in the ranges (vii), (v), and (ii), and these were also the ranges that showed the smallest possible value of determinism (i.e., 0). In terms of the residuals of (iii): they were about 18, as in (ii), which had a null value of determinism, despite the latter being different from 0 for (iii). With the exception of range (iii), all intervals with a determinism other than 0 appeared to show lower residuals. The values of the standard deviation did not seem to be able to explain the goodness of fit. In fact, the Pearson correlation coefficient between the standard deviation and $r_{\mathrm{opt}}$ was $-0.0836063$, whereas between determinism and $r_{\mathrm{opt}}$, it was $-0.6120757$.

7. Conclusions

In this paper, we developed a novel stochastic differential equation (SDE) that fit the model of intraday copper futures prices and was based on the whitening filter approach Wiener (1949), a method characterizing linear time-variant systems. We assumed that the time series of prices could be described by its autocorrelation. This property was obtained through statistical analyses of historical data. Starting from an empirical autocorrelation model obtained from the price time series, we introduced an SDE that was characterized by this autocorrelation and fit the model to the data.

The model we developed produced trajectories close to the real data only in the short run. This was due, on the one hand, to the model fitted by empirical autocorrelation and, on the other, to the predominantly stochastic behavior of the copper price time series, which can only be predicted in the short term Mastroeni and Vellucci (2019, 2022). To investigate the stochastic (vs. deterministic) nature of the data, we followed the intuition behind the recurrence analysis Eckmann et al. (1987); Marwan et al. (2007). In the Appendix A, we also show that the performance of the model improved significantly, when the time series had an empirical autocorrelation function close to the damped cosine model and a determinism measure (DET) based on the recurrence plot close to one. Therefore, future work will include the assumption of more complex functions to model the autocorrelation function of intraday copper futures prices.

Based on the numerical results presented in Section 6.2 regarding the level of determinism through the recurrence plot and the performance of the SDE model, this paper showed that the predictability of copper futures prices is rather weak. In fact, as already said, the model we derived produced trajectories close to the actual data only in the short term. Consequently, the investment risk of copper futures is high, which confirms the results in Mastroeni et al. (2018). Based on the reconstructed copper price phase space, the recurrence plot confirmed that deterministic features coexisted with stochastic ones in copper futures prices. In the intervals considered, the level of determinism was often zero, indicating that the copper futures prices were not predictable.

Compared to other metal futures such as aluminum and nickel, copper futures seem to have a higher investment risk. For instance, the prices considered by Sun et al. (2022) were highly deterministic. In any case, the time series considered here were different in terms of granularity and time span, and intraday copper futures are probably a riskier investment than daily ones. Investigating this aspect could be an interesting issue for a future work.

Additional further research would be to generalize the SDE model to the multivariate case. In fact, since copper prices are correlated with other commodities (such as steel, gold, and oil), it would be interesting to know whether this model could use this information. In any case, co-movements among these series did not arise from the model, but it is necessary to use ad hoc methodologies, e.g., those employed in our previous papers Mastroeni et al. (2021, 2022); Benedetto et al. (2020, 2021), Benedetto et al. (2016). A possible generalization of the SDE model proposed in this paper should integrate these methodologies in the whitening filter approach.