Forecast Relative Error Decomposition

Forecast Relative Error Decomposition

\setstretch{1}

Introduction

Since its introduction in the time series literature [see e.g. Doob (1953), Whittle (1963)], the Forecast Error Variance Decomposition (FEVD) has been largely used in macroeconomics to analyse the effects of shocks and their propogation at different horizons in time [see e.g. Lanne, Nyberg (2016), Isakin, Ngo (2020)]. This approach is easy to implement at the cost of serial and cross-sectional linearity assumptions. It is typically applied on linear dynamic models such as Structural Vector Autoregressions (SVAR) with independent and identically distributed innovations. Moreover, the variance is just a local measure of uncertainty when this uncertainty is small [see Pratt (1963), Arrow (1965) for the definition of risk aversion], and their local interpretation is the basis of the short term mean-variance portfolio management in Markowitz (1953), (2000)]. The goal of our paper is to extend the standard FEVD in the two following directions above: [1] First to account for nonlinear dynamic models (nonlinear serial dependence), that involves asymmetric effects, conditional heteroscedasticity, switching regimes, or extreme risks. [2] Second, to allow not only for the prediction of the future value of the process, but also for the prediction of its nonlinear transforms (nonlinear cross-sectional dependence), in the short, medium and long run. Section 2 provides the background on the variance decomposition formula, its role in defining the standard FEVD, and the local interpretation of this decomposition. In Section 3, we extend the FEVD measure in a nonlinear dynamic framework. We first introduce a general decomposition for the relative forecast errors on a positive transformation of the process of interest, which leads to the so-called Forecast Relative Error Decomposition (FRED). We then apply this decomposition to the updating (i.e. learning) of the forecasts and propose two decompositions based on selected characterizations of the predictive distribution: (i) The Forecast Error Kullback Decomposition (FEKD) is a measure based on the transition density. (ii) The Forecast Error Laplace Decomposition (FELD) is a measure based on the conditional Laplace transform. We demonstrate the suitability of these new decompositions for nonlinear dynamic models in Section 4. We begin with the standard linear (Gaussian) Vector Autoregression (VAR) and compare our measures with the traditional FEVD. Then, we consider different examples in which the nonlinear features are due to the structural interpretation of the process of interest, which can be qualitative, a sequence of count variables, be positive valued, or even in the case of observed volatility-covolatility features be valued as a symmetric semi-positive definite matrix. Numerical illustrations are provided in Section 5, and statistical inference is discussed in Section 6. An application to count observations on cyberattacks is provided in Section 7 and we conclude in Section 8. Technical details are provided in the appendices and online appendices.

Variance Analysis as a Local Approximation

This section introduces two preliminary lemmas that are useful in understanding the Forecast Error Variance Decomposition and its interpretation.

Decomposition of Variance

Let us consider a random vector $Y$ of dimension $n$ and an information set ($\sigma$-algebra) $I$ such that $Z_t$ is measurable with respect to $I_t$.

This decomposition can be written as:

where $'$ denotes a transpose. It relates the multivariate risk of prediction updating (left hand side) to the updating of the variances of the error forecasts (right hand side). The matrix equations (ref)-(ref) can be decomposed by components. They lead to the one dimensional variance decompositions as:

as well as one dimensional covariance decompositions:

in which the information set is not specific to a given component. It is easily checked that these additive decompositions of the variances and covariances do not imply a simple decomposition of the associated correlations $\rho(Y_i,Y_j)=\frac{Cov(Y_i,Y_j)}{\sqrt{\mathbb{V}(Y_i)\mathbb{V}(Y_j)}}$ in terms of the conditional correlations.

Forecast Error Variance Decomposition

We can extend the decomposition of variance to any multivariate stochastic process $(Y_t)$ and its increasing sequence of information sets (or filtration) $(I_t)$, where $I_t=(\underline{Y}_t)$ is the $\sigma$-algebra generated by the present and lagged values of the process. Indeed, the forecast errors at horizon $h$ can be written as a sum of multivariate forecast updates:

By the optimality of the conditional expectations, the forecast updatings are uncorrelated conditional on $I_t$. Then, we get the FEVD.

This decomposition is often written for a strictly stationary process $(Y_t)$, which admits an infinite strong moving average\footnote{Since our objective is an extension of the FEVD to a nonlinear dynamic framework, it is important to consider strong moving average models, that are defined from strong white noise instead of weak moving average models where the noise is just assumed to be zero mean, with fixed variance and no serial correlation.} representation of the form:

where $\varepsilon_t$ is a strong white noise, that is a sequence of i.i.d random vectors, with zero mean and variance-covariance $\Sigma$, and the moving average coefficients satisfy the square integrability condition $\sum_{j=0}^{\infty} \vert\vert A_j \vert \vert^2 < \infty$. When this moving average representation is invertible, the process also admits an (infinite) autoregressive representation:

and the information generated by the current and lagged values of the process $(Y_t)$ is equal to the information generated by the current and lagged values of the noise $(\varepsilon_t)$, that is:

Then, the FEVD becomes:

since:

In this special moving average case with iid noise, the decomposition is path independent, that is independent on the values of process $(Y_t)$ before time $t$. These decompositions in the strong moving average case, that is in the strong linear dynamic case, have been initially analyzed in a systematic way by Whittle (1963) [see also Doob (1953)].

Local Analysis of Risk

The comparison of multivariate quantitative risks (i.e. the losses) $X$ and $Y$ is based on the notion of stochastic dominance at order 2. $Y$ is riskier than $X$ if and only if $\mathbb{E}[\upsilon(Y)] \geq \mathbb{E}[\upsilon(X)]$, for any increasing convex function\footnote{This stochastic dominance can be equivalently defined in terms of preference, with now utility functions $u$ that are increasing concave.} $\upsilon$ (such that the expectations exist) [see Rothschild, Stiglitz (1970), Vickson (1975), Fishburn, Vickson (1978).]. In the univariate case, it is well known that the quadratic function $y \rightarrow y^2$ that underlies the definition of variance is convex, but not increasing (except if $X$ and $Y$ are nonnegative). Then the variance is not directly appropriate for measuring risk. However, we have the following local expansion.

Proof: (i) This is a consequence of Jensen's inequality. \\

(ii) Indeed we have:

and by taking the expectation of both sides:

Since $Y$ is assumed close to zero, the variance-covariance matrix is seen as a local risk measure for small risks and by construction does not account for asymmetric risks (which would require an expansion up to order 3), or to extreme risks (which would require an expansion up to order 4). This second--order expansion has been the basis for defining the local version of absolute risk aversion\footnote{In the multidimensional framework, the risk aversion is a matrix directly linked to the Hessian at 0 of function $v$, that is, $\frac{\partial^2v(0)}{\partial y \partial y'}$.} [Arrow (1965)], or for justifying the mean-variance portfolio management in finance [Markowitz (1952), (2000)]\footnote{In the standard financial application with $Y$ as a vector of asset returns, the interest is in the risk of the portfolio returns $\alpha'Y$, where $\alpha$ is the vector of portfolio allocations in values. Then the risk becomes scalar, that is we have $\upsilon(Y)=u(\alpha'Y)$, where $u$ is defined on $\mathbb{R}$. We deduce that the matrix of risk aversion $\frac{\partial^2\upsilon(0)}{\partial y \partial y'}=\frac{d^2u(0)}{dw^2}\alpha \alpha'$ is of rank 1. It involves the effect of portfolio allocation and a scalar risk aversion.}. Since the risk on a price evolution $y=p_{t+h}-p_t$ increases generally with the term, the mean-variance management is appropriate in the short run, that is with frequent predictive updating, and frequent learning.

Forecast Relative Error Decomposition (FRED)

As mentioned in section 2.3, the standard FEVD has two drawbacks: (i) While it is appropriate for a local analysis of risk, it is not appropriate for asymmetric risks, or extreme risks. (ii) Moreover it is usually applied to pointwise predictions of $Y$, not to predictions of nonlinear transformations of $Y$. The aim of this section is to provide decompositions that are more appropriate for prediction and learning in a nonlinear dynamic framework. We first provide a general FRED decomposition. This decomposition is applied to one dimensional positive transformations of the process, that are either transition densities, or conditional Laplace transforms of process $(Y_t)$, leading to the so-called Forecast Error Kullback Decomposition (FEKD) and Forecast Error Laplace Decomposition (FELD), respectively.

Relative Forecast Updating

Let us consider a univariate positive process $(Z_t)$ and the sequence $(I_t)$ of increasing information sets such that $Z_t$ is measurable with respect to $I_t$. We can construct an analogue of the FEVD by considering the relative forecast update:

Then we have:

By taking the log of both sides and then the conditional expectation given $I_t$, we get the Forecast Relative Error Decomposition (FRED):

Proof: The nonegativity is a consequence of Jensen's inequality. For instance, let us consider the left hand side. Then:

The proof is similar for the terms in the right hand side, noting that $\mathbb{E}(\cdot | I_t) = \mathbb{E}(\mathbb{E}(\cdot |I_{t+h})|I_t)$ by the Law of Iterated Expectation. \qed \\

The FRED is locally a variance decomposition written on the relative error forecasts. Indeed we have locally:

and

since $\mathbb{E}\left[\frac{Z_{t+h}-\mathbb{E}(Z_{t+h}|I_t)}{\mathbb{E}(Z_{t+h}|I_t)}\bigg\vert I_t\right]=0$. Therefore, the main difference between the FRED and the FEVD is that locally the variance expansion is written on the relative forecast errors instead of the absolute forecast errors and can be applied to any positive transformation $Z_t$ of the multivariate process $(Y_t)$. The FRED (ref) can be written as:

where the left hand side of equation (ref) $\gamma(h|I_t)$, say, measures the risk on the prediction errors at horizon $h$, and the generic term of the right hand side, $\gamma(k,h|I_t)$, has a forward interpretation. More precisely, we have:

Therefore, $\gamma(k,h|I_t)$ is the expectation at date $t$ of the risk on the short run prediction error of $Z_{t+h}$ at date $t+k$. As usual, this forward interpretation involves three dates: $t$, $t+k$ and $t+h$.

Forecast Error Kullback Decomposition (FEKD)

To account for cross-sectional nonlinearities and avoid pointwise predictions, a first solution is to consider the transition predictive densities. Let us assume a strictly stationary Markov process $(Y_t)$ with transition density at horizon $h$ and time $t$ given by $f(y,h|Y_t)$. These densities are positive transforms of $Y_t$. They are well defined for horizons $h\geq 1$, but not at horizon $h=0$, where the value $Y_t$ is perfectly known and the predictive distribution degenerates into a point mass at $Y_t$. Up to this degeneracy, we deduce from Lemma 4 the Forecast Error Kullback Decomposition below:

We get a decomposition that involves measures of risk on the updating of predictive densities. When $Y$ is a continuous variable, the FEKD is invariant by one-to-one differentiable transformations of $Y$. Indeed, due to the ratios, the Jacobian effect disappears in the decomposition formula. The generic term on the right hand side of the FEKD is:

that is the conditional Kullback proximity measure (or contrast, or divergence) between the two conditional densities. \\

A main difference between the FEKD and FEVD is that the decomposition (ref) can be written for any value of $y$, whereas there is a single FEVD. Typically, let us consider a univariate stationary process $(Y_t)$. From $T$ observations $Y_1,...,Y_T$, we can derive the sample distributions, then the sample deciles, and compare the FEKD evaluate at difference deciles. This will provide a more detailed analysis of risk with clear applications to the dynamic analysis of inequality, or the analysis of financial risks (where these quantiles are usually called Value-at-Risk (VaR)). \\

When the process $(Y_t)$ has a dimension larger than 2, the notion of quantile has not yet been defined, and it is less clear how to select an appropriate grid of values $y$ to which the FEKD will be applied. A solution can be to apply the approach by deciles to combinations of components of $Y_t$, if such combinations have an economic meaning. Combinations are more easily treated by means of Laplace transforms, as seen in the next subsection. Note that the FEKD can exist in cases where the FEVD has no meaning. Indeed, the FEVD requires that the observed $Y_t$ are square integrable, so it cannot be applied to data with fat tails. We give below in Section 4.1.2 the example of the Cauchy AR(1) model, where the FEKD exists, but not the FEVD.

Forecast Error Laplace Decomposition (FELD)

Let us consider a process $(Y_t)$ of dimension $n$. Its conditional Laplace transform at horizon $h$ is defined by:

where $u\in D \subset \mathbb{R}^n$ and such that the expectation exists on domain $D$. It is known that the knowledge of the Laplace transform is equivalent to the knowledge of the conditional distribution for the Gaussian case, or if the process satisfies some positivity restrictions [see Feller (1971) and the examples in Section 4]. Then we can apply the FRED to these nonlinear transformations.

As for the FEKD, the FELD (ref) defines several decompositions: $\gamma(h|u,I_t)=\sum_{k=0}^{h-1}\gamma(k,h|u,I_t)$, as many as selected combinations $u$, $u\in D$. In some applications such combinations can have economic interpretations to define portfolio allocations when $Y_t$ is a vector of $n$ asset prices at date $t$, or to combine income and wealth in the analysis of inequality. When $(Y_t)$ satisfies some “positivity" restrictions, the Laplace transform with positive argument $u$ characterizes the distribution. Moreover, we have $\mathbb{E}\left[\exp(-u'Y_{t+h})\vert I_t\right]\leq 1$, and the FELD always exists, even if the distribution of $Y_t$ has very fat tails.

Decomposition of Risk Premium

The FELD can be interpreted as a decomposition of risk premiums for spot prices. Consider a decision maker with exponential utility function $u(y)=-\exp(-uy)$. The certainty equivalent $\pi(u)$ is the level of wealth which makes the decision maker indifferent from the expected outcome of a lottery on $Y$. That is:

which is a function of parameter $u$, which is the Arrow-Pratt scalar measure of risk aversion. Equivalently, for a risky asset with price $Y_t > 0$, we can write:

which is the spot value of the asset at time $t$ and horizon $h$. Intuitively, it is a contract written at time $t$, which captures the value of delivering the asset at some future horizon $h$.\\

For an investor with risk aversion parameter $u$, the FELD in (ref) can be written as:

Let us now consider the economic interpretations of this decomposition formula. The term $\pi(u,h|I_t) - \mathbb{E}(Y_{t+h}|I_t)$ is the difference between the value (price) of $Y_{t+h}$ at date $t$ and its historical conditional expectation. This difference is positive (by Jensen's inequality) and usually interpreted as a risk premium. Therefore, (ref) provides a decomposition of this risk premium. More precisely, the term $\pi(u,k,h|I_t)=\mathbb{E}\left[\pi(u,h-k|I_{t+k})\bigg \vert I_t \right]$ is the value of a forward contract of the asset, written at time $t$, for a payment at time $t+k$ and delivery at time $t+h$. Hence, the generic term on the RHS of (ref) captures the difference in values $\pi_f$ of the forward contracts for payment at time $t+k$ and $t+k+1$ for the delivery of the asset at time $t+h$. As $h$ varies, we get a decomposition of the term structure for the risk premium, or equivalently of the spot value (price) as:

that is compatible with the no dynamic arbitrage condition between spot and forward contracts. There is a debate on the valuation approach to be chosen for contingent assets [see e.g. Embrechts (2000)]. For financial assets traded on very liquid markets, this is usually done by introducing a stochastic discount factor to satisfy the no dynamic arbitrage opportunity condition. The situation is different for individual insurance contracts or for operational risks [see the discussion of cyber risk in Section 7]. The certainty equivalent principle is a more appropriate valuation approach in such frameworks and we have checked ex-post that it is compatible with the no arbitrage opportunity assumption.

Examples

In this section, we consider different dynamic models for which we derive closed form decompositions for either the FEKD, or the FELD.

Examples of FEKD

The Gaussian VAR(1)

Let us assume that the $n$-dimensional stationary process $(Y_t)$ satisfies:

where the eigenvalues of $\Phi$ have a modulus strictly smaller than 1 and the $\varepsilon_t$'s are i.i.d. Gaussian $\varepsilon_t \sim N(0,\Sigma)$. Then, the conditional distribution of $Y_{t+h}$ given $Y_t$ is Gaussian with mean $\Phi^hY_t$ and variance-covariance matrix $\Sigma_h = \Sigma + \Phi \Sigma \Phi'+...+\Phi^{h-1}\Sigma(\Phi')^{h-1}$. The following proposition provides the closed form FEKD for the Gaussian VAR(1).

Proof: See the Appendix A.1.1. \\

For the Gaussian VAR(1) the FEVD is:

The FEKD in Proposition 3 differs from the FEVD, since it depends on $Y_t$ and on the argument $y$. Its functional quadratic form in $y$ provides information on the prediction errors and their decomposition in the tails when we focus on the components in $y$ and the squares. The decomposition of the intercept is just to balance the change in tails, since all the predictive distributions are unit mass. The component in $y$ depends on $Y_t$, which corresponds to the forward interpretation of the elements in the decomposition. The quadratic element is independent of $Y_t$ due to the conditional homosecedasticity of the Gaussian VAR(1) model. This element is written on the inverses of the prediction variance, not on the prediction variances as in the FEVD. This difference is the analogue of the two equivalent filters, introduced by Kalman for the linear state space model. Indeed, the standard covariance filter is based on the direct updating of $\Sigma_h$, whereas the information form of the filter coincides with the direct updating of the inverse $\Sigma^{-1}_h$ (called information).

The Cauchy AR(1)

Let us consider the stationary univariate process $(Y_t)$ defined by:

where $(\varepsilon_t)$ is a Cauchy distributed strong white noise. This distribution admits the density $f(\varepsilon)=\frac{1}{\pi}\frac{1}{1+\varepsilon^2}$, and its characteristic function is given by: $\mathbb{E}(\exp(iu\varepsilon))= \exp(-|u|)$, where $i$ is the imaginary number. Then we have:

where:

Therefore the conditional distribution of $Y_{t+h}$ given $Y_t$ is a Cauchy distribution with drift $\varphi^hY_t$ and scale $\sigma \frac{1-|\varphi|^h}{1-|\varphi|}$. We deduce that:

and then:

This is the generic term of the RHS of the FEKD. In this framework the conditional variance at horizon $h$, that is $\mathbb{V}\left(Y_{t+h}|Y_t\right)$, does not exist and the FEVD does not exist as well. However, the elements $\gamma(k,h|I_t)$ in the decomposition above exist, since the transformed variable $\log(\alpha+\beta\varepsilon_t^2)$ is integrable with respect to the Cauchy distribution.

Markov Chain

Let us start by considering a stationary Markov chain with two states $0,1$, or equivalently, a stationary Markov binary time series $(Y_t)$. The transition matrix of the chain can be parameterized by the marginal probability $\pi = P(Y_t =1)$ and a persistence parameter $\lambda$, $0 \leq \lambda < 1$. We have:

Then by iterated expectations we deduce:

Proof: See Appendix A.1.2. \\

Such univariate binary processes have attracted considerable attention to describe the business cycle recession/expansion periods [see e.g. Estrella, Mishkin (1998), Chauvet, Potter (2005), or Kauppi, Saikkonen (2008)]. \\

This term can be compared with the generic term in the standard FEVD.

Proof: See Appendix A.1.3. \\

Now we can deduce the extension to the general case. Suppose $Y_t$ is a stationary Markov chain with $n$ possible states. Let $X_t = (X_{1,t},...,X_{n,t})'$, with $X_{i,t}=1$, if $Y_t$ is in state $i$, and zero otherwise, for $i=1,...,n$. The knowledge of $Y_t$ is equivalent to the knowledge of $X_t$, which can take on the values: $(1,0,...0)'$, $(0,1,0,...,0),...,(0,...,0,1)'$. Let us characterize its $h$-step transition probabilities with the matrix $P^{h}$, whose elements are defined as $p^{(h)}_{ij} = P(Y_{t+h}=i|Y_t=j)$ for all $t\geq0$. Then, the FEKD is given by the following proposition.

Proof: See Appendix A.1.4. \\

The result presented in Propostion 4 is a special case of Proposition 5.

Proof: See Appendix A.1.5. \\

Examples of FELD

The FELD has a simple form and is easy to interpret for the dynamic affine model, which is called the Compound Autoregressive (CaR) model in discrete time [see Duffie, Filipovic, Schachermayer (2003), Darolles, Gourieroux, Jasiak (2006)]. We first review the CaR models and some of their dynamic properties. Then, we discuss in detail the application to Gaussian processes, Integer Autoregressive (INAR) models, Negative Binomial Autoregressive (NBAR) models, Markov Chains, Autoregressive Gamma and Wishart processes. Remember that the FELD has an interpretation as a decomposotion of risk premium where $Y$ is a value (see Section 3.4).

Dynamic Affine Model

The process is assumed Markov of order 1 with a conditional log-Laplace transform which is affine in the conditioning value\footnote{The extension to a Markov of order $p$ is straightforward.} $Y_t$. Thus, the Laplace transform at horizon 1 can be written as:

where $c(u) =\log \mathbb{E}[\exp(-u'Y_t)]$ is the unconditional log-Laplace transform of $Y_t$ and the function $a(\cdot)$ captures all nonlinear serial dependence features. The affine property remains satisfied at any forecast horizon. More precisely, we have:

where $a^{\circ h}(\cdot)$ is function $a(\cdot)$ compounded $h$ times with itself. Then, the FELD becomes:

This decomposition involves the conditional expectations of $Y_{t+k}$ given $Y_t$ and can be written as:

Whereas the FELD involves the conditional expectations only, the nonlinear dynamic features are taken into account by the “weightings" of the expectation that depend on the function $a(\cdot)$, that summarizes the nonlinear serial dependence. \\

It is known that the log-Laplace transform admits a Taylor expansion in terms of cumulants in a neighborhood of $u=0$. In particular, the first-order expansion of the conditional log-Laplace transform shows that the conditional expectation is affine in the conditioning variable $Y_t$. More precisely, we have:

and

Then we deduce a closed form FELD for dynamic affine models.

Proof: See Appendix A.2.1. \\

Therefore, we get a decomposition of the type:

or equivalently decompositions of the functions $\alpha(h,u)$ and $\beta(h,u)$ into $\sum_{k=0}^{h-1}\alpha(h,k,u)$ and $\sum_{k=0}^{h-1}\beta(h,k,u) $ respectively. These decompositions depend on both the term $h$ and the argument $u$. The decomposition of function $\alpha(h,u)$ is especially appealing due to its interpretation in terms of nonlinear Impulse Response Functions (IRF). Indeed, let us consider a given shock of magnitude $\delta$ on the level $Y_t$. Note that in our nonlinear framework, the magnitude of the shock $\delta$ is constrained by the domain of $Y_t$. It can be multivariate if $Y_t$ is multivariate, constrained to be either 0, or -1 (resp. 0, or 1) if $Y_t$ is binary with value 1 (resp. value 0), and so on. We do not discuss the sources of the shock and if they are identifiable or controllable. The effect on the predictive distribution at horizon $h$ is $\alpha(h,u)'\delta$. Compared to the standard linear approach of the IRF, we see that the IRF depends on the argument $u$. In other words, this measure changes with the preference (risk aversion) of the analyst. Moreover, for a given $u$, the decomposition will change, since in a nonlinear dynamic forecast, the spot and forward short run updatings of the predictive distributions will differ. \\

The FELD can be easily compared to the FEVD. Let us consider the one-dimensional case for expository purposes. In an affine dynamic model, the conditional variance $\mathbb{V}(Y_{t+h}|I_{t+k})$ is an affine function of $Y_{t+k}$, and thus the generic term in the FEVD (ref) is also an affine function of $Y_t$. More precisely, we get:

Proof: See Appendix A.2.2.

Strong Linear VAR(1) Model

Let us consider the strong VAR(1) model:

where the $\varepsilon_t$'s are i.i.d with the log-Laplace transform:

The conditional Laplace transform is:

When the eigenvalues of $\Phi$ have a modulus strictly smaller than 1, we can write the infinite moving average representation of process $(Y_t)$ as:

Therefore, the unconditional log-Laplace transform of process $(Y_t)$ is:

The strong VAR(1) is an affine model, with $a(u)=\Phi'u$ and $a^{\circ h}(u)=(\Phi')^hu$. Since the dynamic is linear and the function $a(\cdot)$ is also linear, the FELD is greatly simplified. We get:

Due to the linear dynamic, the FELD does not depend on the value $Y_t$ of the conditioning variable. However, there are still effects on the decomposition of the cross-sectional heterogeneity, that is the non-Gaussian distribution of the errors $(\varepsilon_t)$. If the error is Gaussian $\varepsilon_t \sim IIN(0,\Sigma)$, we get $b(u)=\frac{u'\Sigma u}{2}$. Then, the right hand side of the decomposition becomes:

We recover the right hand side $\sum_{k=0}^{h-1} \Phi^k \Sigma \left(\Phi'\right)^k$ in the FEVD for the VAR(1) model (see equation (2.8)). To summarize:

Markov Chain

Let us return to the Markov chain example discussed in section 4.1.2 with $n$ states. The conditional log-Laplace transform is given by:

where $u=(u_1,...,u_n)'$ and $\widetilde{\log}(A)$ is a matrix whose elements are the logged elements of matrix $A$. Then the FELD is given by the following proposition:

Proof: See Appendix A.2.3.

INAR Model

For expository purposes, let us consider the Integer Autoregressive model (INAR) of order 1 introduced by McKenzie (1985), Al-Osh, Azaid (1987). The process is defined by:

with $\mathcal{B}_t(p) \circ Y_{t-1}=\sum_{j=1}^{Y_{t-1}}U_{j,t}$, where the variables $\varepsilon_t$, $U_{j,t}$, $j,t$ varying, are independent, $U_{j,t}$ follows the same Bernoulli distribution $\mathcal{B}(1,p)$ and $\varepsilon_t$ the Poisson distribution $P(\lambda)$, with $0 \leq p < 1$ and $\lambda\geq0$. By convention, $\sum_{j=1}^{0}U_{j,t} = 0$. The INAR model is a CaR model, with $a(u) = -\log\left[p\exp(-u)+1-p\right]$, $c(u)=\frac{-\lambda}{1-p}\left[1-\exp(-u)\right]$. In particular, the marginal distribution of $Y_t$ is Poisson $P\left(\frac{\lambda}{1-p}\right)$. It is easily checked that: $a^{\circ h}(u)=-\log\left[p^h\exp(-u)+1-p^h\right]$.

Proof: See Appendix A.2.4. \\

Cox, Ingersoll, Ross and Autoregressive Gamma Process

Let us consider the univariate Autoregressive Gamma Process of order 1, which is the time discretized Cox, Ingersoll, Ross model [see Gourieroux, Jasiak (2006)]. This dynamic model is the benchmark for the dynamic analysis of short run interest rates, or for one dimensional stochastic volatility. This is an affine model where the transition distribution is constructed from a gamma distribution with path dependent stochastic degree of freedom. It depends on two parameters: $\beta\geq 0$, $\delta\geq0$, and the Laplace transform corresponding to $a(u)=\frac{\beta u}{1+u}$, $c(u)=-\delta \log\left[1+\frac{u}{1-\beta}\right]$. In particular, the marginal distribution of the process is such that: $(1-\beta)Y_{t} \sim \gamma(\delta)$. It is easily checked that $a^{\circ h}(u)=\frac{\beta^h u}{\left[1+\frac{1-\beta^h}{1-\beta}u\right]}$ and that $\frac{da(u)}{du} = \beta$. Therefore, we get the following FELD decomposition of the functional IRF (the decomposition of the intercept $\beta(h.u)$ is provided in Appendix A.2.5).

Proof: See Appendix A.2.5. \\

Hence, the proportion $\frac{\alpha(h,k,u)}{\alpha(h,u)}$ is given by:

As $u$ increases, the numerator is a positive decreasing function, while the denominator is an increasing function of $u$. Hence, the proportion $\frac{\alpha(h,k,u)}{\alpha(h,u)}$ is a decreasing function of $u$ as well.

The Wishart Process

The analysis of risk measured by volatility can be extended to the multivariate framework and volatility-covolatility matrices. This leads to the Wishart Autoregressive (WAR(1)) process that is the multivariate analogue of the ARG(1) process [see Gourieroux, Jasiak, Sufana (2009) for discrete time and Cuchiero et al. (2011) for continuous time]. Due to the matrix framework and the positivity satisfied by a volatility-covolatility matrix, the conditional Laplace transform, is usually written as:

where $\Gamma$ is a matrix of arguments assumed symmetric positive semi-definite and $\operatorname{Tr}$ denotes the trace operator that sums up the diagonal elements of a square matrix. Then we get:

Proof: See Appendix A.2.6. \\

It can be checked that this element reduces to the element in the FELD of the ARG (see Section 4.2.5) in the one dimensional case.

Illustrations

We now provide some empirical illustrations of the theory presented in the preceding section.

FEKD for Gaussian VAR(1) Process

We consider a bivariate VAR(1), $Y_t=(Y_{1,t},Y_{2,t})'$, with autoregressive parameter $\Phi =

$. The eigenvalues of matrix $\Phi$ are $\lambda_1=0.7$ and $\lambda_2=0.4$, so the model admits a stationary solution for $(Y_t)$. Let the covariance matrix of $\varepsilon_t$ be given by $\Sigma =

$, that is, the innovations have unit variance and are positively correlated. \\

The first illustration of interest is how the total FEKD, that is, the left hand side of (ref), varies with different values of $(Y_t,y)$. In fact, this total relies on how “far" the value $y$ is from the mean of the conditional predictive density at time $t$, which is equal to $\Phi^hY_t$. To see this, let us define the Mahalanobis distance between $y$ and $\Phi^hY_t$, given by:

which is an extension of the standard deviation in a multivariate setting. For expository purposes, suppose $Y_t = (2,1)'$ and $h=10$. Then, $\Phi^h(2,1)'=(0.028,0.56)'$. We plot in Figure 1 the relationship between the total FEKD in (ref) and the distance $d(y,\Phi^hY_t)$. The total FEKD increases exponentially with respect to $d(y,\Phi^hY_t)$. This means that values of $y$ further from the mean $\Phi^hY_t$ will have a larger marginal increase in total FEKD. Moreover, note that two points with the same $d(y,\Phi^hY_t)$ will also have the same total FEKD. For example, the points $y=\Phi^h(2,2)'=[0.038, 0.075]'$ and $y=\Phi^h(2,0)'=[0.019, 0.038]'$ both have a Mahalanobis distance of 0.0101 and a total FEKD of 0.1415. However, even if two points share the same Mahalanobis distance, their decompositions can be very different. Let us return to the values $y=\Phi^h(2,2)'=[0.038, 0.075]'$ and $y=\Phi^h(2,0)'=[0.019, 0.038]'$. \\

In the first two graphs of Figure 2, we illustrate the decompositions of $b(h|Y_t)y$ and $y'c(h|Y_t)y$ respectively. It is clear that the point $\Phi^h(2,2)$ has stronger effects on both components. However, when we consider the total sum $b'(h|Y_t)y+y'c(h|Y_t)y$, it is not surprising that both $y=\Phi^h(2,2)'=[0.038, 0.075]'$ and $y=\Phi^h(2,0)'=[0.019, 0.038]'$ produce an identical number, as shown in the bottom left graph. Indeed, since the total FEKD is the same for both points and component $\alpha(h|Y_t)$ is independent of $y$, the sum $b'(h|Y_t)y+y'c(h|Y_t)y$ must also be the same for two points with the same Mahalanobis distance. Furthermore, we note that the comparison of decomposition for these two points is somewhat inconsequential in understanding the dynamics of the VAR(1), since in the bottom right graph, we see that the sum $b'(h,k|Y_t)y+y'c(h,k|Y_t)y$ contributes a very small fraction of the total FEKD in comparison to $a(h,k|Y_t)$. \\

On the other hand, this latter remark is not always valid for all $y$. Indeed, let us now illustrate an example of points that have a much higher Mahalanobis distance from the mean $\Phi^hY_t$, that is, in the tail ends of the predictive density. For instance, consider $y=\Phi^h(2,-99)'=[-0.91, -1.83]'$ and $y=\Phi^h(2,101)'=[0.97, 1.94]'$, which have a Mahalanobis distance of 1.0146. \\

In Figure 3, we can see again that the decomposition of components $b'(h|Y_t)y$ and $y' c(h|Y_t)y$ can be very different, even if two points share the same Mahalanobis distance. Also, note that the components of $y' c(h|Y_t)y$ are much larger in magnitude than the components of $b(h|Y_t)y$, which suggests that $y' c(h|Y_t)y$ (the quadratic term) is much more informative on dynamic behaviour in the tails. Furthermore, as seen in the last graph of Figure 3, the sum $b'(h|Y_t)y+y'c(h|Y_t)y$ now takes up a much larger fraction of the total FEKD. Hence, the decomposition is much more sensitive to the choice of $y$ in the tails than it is near the mean. \\

FELD for INAR(1) Model

Let us consider the INAR(1) process and its FELD derived in Corollary 6, eq. (ref). The decomposition depends on the current value $Y_t$, on the mean of the Poisson innovation $\lambda$, on the persistence parameter $p$, on the risk aversion $u$, and on the horizon $h$. For expository purposes, we set $Y_t=3$ and $\lambda=2$ in the illustrations, noting that they only enter the FELD as scaling constants and do not change the properties of the decomposition. \\

In Figure 4, we present graphs depicting the relationship between the total FELD and the risk aversion parameter $u$ in $[0.1,2.8]$, under different values of persistence parameter $p$. The lines in each figure represent the FELD for different horizons, with darker shades representing later horizons. For instance, the black line represents the total FELD for horizon 10. \\

In general, the total FELD is an increasing function of $u$. That is, when a decision maker is more risk averse, the total information gain at a future horizon is also higher, measured by means of the relative change in the Laplace transform of the process. The FELD is also influenced by the level of persistence in the INAR(1). When persistence is low, the information gain in each horizon is similar; in the bottom right graph for $p=0.05$, the line for each horizon is almost stacked on top of each other, whereas the lines in the later horizons are higher in the top left graph for $p=0.99$. Indeed, persistence also influence the steepness of the FELD in each horizon. \\

In Figure 5, we consider decompositions of the FELD at horizons $h=1,...,10$, based on a grid of values in $u=(0.5,1,2)$ and $p=(0.1,0.5,0.95)$. The height of each bar in the graphs correspond to the value of the FELD on the left-hand side of (ref). They are decomposed into smaller bars with different shades, with the darker shades corresponding to a larger $k$ on the right-hand side of (ref). \\

The risk aversion parameter only influences the total FELD amount, but it has no impact on the decomposition of its components. Instead, the persistence of the process influnces the FELD in two ways. Firstly, at higher levels of persistence, there is higher contribution of short term updates in previous horizons. For example, let us consider the FELD for horizon 10 at $u=2$ (i.e. the graphs in the third row). When there is almost no persistence ($p=0.1$), the bar at horizon 10 is completely in black, meaning that only the update from horizon 9 to 10 contributes to the total FELD. When persistence is high ($p=0.95$), the bar at horizon 10 contains various shades of grey, which means that updates from previous horizons (e.g. 1 to 2, 2 to 3, 3 to 4, etc.) also contribute to the total FELD. This is also depicted in Figure 6, where the FELD is instead taken as a fraction out of 100. \\

Secondly, for lower levels of persistence, the total FELD seems to “converge" towards a long-run value. This is due to the stationarity of the INAR(1) process. Indeed, although it is not visible in this illustration, the FELD also converges to a long-run value for $p=0.95$, albeit at a much further horizon. We plot in Table 1 the limiting values of the FELD as $h$ tends to infinity for various combinations of $(p,u)$. To no surprise, the limiting values are increasing in both $u$ and $p$. \\

Another decomposition of interest is the separation of the FELD into the time-varying component $\alpha$ and the fixed component $\beta$. This is exemplified in Figure 7 below for the INAR(1) process with $(p,u)=(0.7,3)$. Although the FELD is increasing over horizon $h$, it is the fixed component that plays a gradually larger role in comparison to the time varying component. Indeed, we expect to see this behaviour of the $\alpha$ component, since it has the interpretation of an impulse response function; for stationary processes, the IRF should exhibit transitory behaviour rather than long-run or permanent changes. Looking at a more granular decomposition of the $\alpha$ and $\beta$ componenents themselves, there are no special patterns that differ greatly from the overall FELD. \\

FELD for ARG(1) Process

Let us now illustrate the properties of FELD under the ARG(1) process, with focus on the terms $\alpha(h,u)$ and $\alpha(h,k,u)$. In Figure 8 below, we present the relationship between the total FELD term $\alpha(h,u)$ and risk aversion $u$ for varying levels of parameter $\beta$, which characterizes the serial dependence. \\

Let us first consider the cases where $\beta$ is relatively high, say above $\beta = 0.9$. For all horizons $h$, the total FELD term $\alpha(h,u)$ is an increasing function of $u$. However, the rate of increase differs across horizons. In particular, when $u$ is small, the total FELD is higher for the later horizons $h$ (i.e the darker lines are above the lighter ones). This relationship reverses when $u$ is high and there seems to be a point at which the lighter lines “cross" the darker ones. When $\beta$ is low, this relationship is seemingly absent, and the lighter lines are always higher than the darker ones. However, there is still a “cross" point for these graphs, but due to the scaling of the axis, they cannot be seen clearly, since they appear at very low values of $u$. Indeed, the value of $\beta$ influences the crossing condition between horizons. In particular, it can be shown that the crossing condition for horizons $h$ and $h+1$ is given by:

In the left graph of Figure 9, we depict the curves for $u^*(h,\beta)$ as a function of $\beta$. Each line corresponds to an horizon $h$, with the darker lines representing later horizons. We can see that for smaller $h$, the total FELD will cross with $h+1$ at a much higher level of risk aversion $u$. For instance, at $\beta=0.90$, the curve for $h=1$ will cross with $h=2$ at $u*=3.737$. However, for the curve $h=2$ and $h=3$, the crossing point is at $u^*=1.047$. These values can be seen on the right graph of Figure 9. \\

Finally, we consider the properties of the FELD decomposition for the ARG(1) process. We consider a grid of values on $(u,\beta)$, for $u=0.5,1,3$ and $\beta=0.1,0.5,0.9$ in Figure 10, which shows the fraction $\frac{\alpha(h,k,u)}{\alpha(h,u)}$. There are two main takeaways from this demonstration. Firstly, as the autoregressive parameter $\beta$ increases (i.e going from top to bottom), then the contribution of previous horizons to the current horizon is higher. Secondly, as the risk aversion parameter $u$ increases (i.e. going from left to right), the contribution of previous horizons to the current horizon is lower. Hence, at the bottom left graph, we can see that even for horizon $h=10$, there is significant contribution from $\alpha(h,k,u)$ for values $k=1,...9$. On the other hand, at the top right graph, the bar for horizon $h=10$ is almost completely homogenous in colour. \\

Statistical Inference

Let us focus on statistical inference for the FELD, the approach being similar for the FEKD. The decomposition is functional, indexed by the current value $Y_t =y$ and the risk aversion parameter $u$. Then we consider the functional estimator in both the parametric and nonparametric frameworks.

Parametric Dynamic Models

When the dynamic model is parameterized as in the examples of Section 4.2, the decomposition takes the form:

where $\theta$ denotes the parameter, $u$ captures the risk aversion and $y$ is the generic value of $Y_t$. If $\hat{\theta}_T$ is a consistent estimator of $\theta$ and is asymptotically normal:

the FRED can be estimated by plugging in $\hat{\theta}_T$ in the decomposition:

Moreover, uniform confidence bands can be derived by applying the $\delta$-method with respect to $\hat{\theta}_T$. Indeed, the doubly indexed vector $\text{vec}\ \hat{\gamma}_T(u,y)=\left[\hat{\gamma}_T(k,h,\vert u, y), \ h=1,...,H, \ k=0,...,h-1\right]$ is asymptotically normal. \\

Nonparametric Dynamic Model

Nonparametric inference can also be applied under the Markov assumption in the one-dimensional case. Let us for instance consider the FELD. The estimation approach is in four steps:

Application to Cyberrisk

The prevalance of the internet in our daily lives means that individuals are now able to communicate, transfer and store large amounts of information with just a mobile device. This has significantly improved the way in which businesses operate on a daily basis and has transformed the structure of our modern economy. For example, hospitals have adopted digital patient records which can be accessed by any institution in their network, and alleviated the need to maintain or deliver physical patient files. However, this also offers an opportunity for bad faith actors to intercept or steal information, leading to potentially disastrous outcomes for the victims involved. As such, there is demand for businesses and government agencies to model and quantify the risk of cyber attacks in order to insure against these prospects. In this section, we demonstrate how the FELD can be used in the multivariate Negative Binomial Autoregressive framework in studying the decomposition of frequency in cyber attacks.

The Challenges and the Available Data

The data on cyber attacks is difficult to obtain. Firstly, there is no consensus on what a cyber attack is. Indeed, there are many definitions, and in a sense it is an umbrella term which includes a variety of illicit behaviours or acts to obtain digital information. Secondly, the collection of data on cyber attacks is rather limited. Although the internet is available almost everywhere today, its accessibility was much less so even just 20 years ago. Hence, tracking cyber attacks has eased only in recent years. Moreover, recording reliable and confirmed cyber attacks is a challenging task, since firms or organizations that are subject to these attacks have an incentive to hide their occurrences. Nonetheless, recent research has appealed to the Privacy Rights Clearinghouse (PRC) dataset, which includes information on publicly reported data breaches across the United States between 2005 and 2022\footnote{https://privacyrights.org/data-breaches. Other cyber databases are Advisen and SAS Oprisk [see Eling, Ibragimov and Ning (2023) for a comparison]}. The PRC was founded in 1992 by the University of San Diego School of Law. The data are gathered from different sources including the Attorney General offices, government agencies, nonprofit websites and media. The reports do not follow a consistent procedure, which may lead to a lack of accuracy and representativeness of the data. Nevertheless, this database is usually employed to analyze cyberrisk [Eling and Loperfido (2017), Eling and Jung (2018), Barati and Yankson (2022), Lu et al. (2024)]. \\

We focus our attention on the modelling of cyber attack frequency counts and adopt the sample used in Lu et al. (2024). The data contains four types of breaches defined as by the PRC:

The plots of these four time series are shown in Figure 12 below.

All four series feature a large number of zero counts, with 146, 138, 440 and 108 instances for DISC, HACK, INSD and ELET respectively. They also have varying levels of occurrences; for instance, while HACK and INSD can both be considered breaches with malicious intent, the latter occurs much less frequently on average. This is due to the fact that it is much easier for an outsider to gain access without getting caught (such as using an IP spoofer), than an insider to attempt a breach. To gain insights on the distribution of each series we provide summary statistics and density plots of each type of breach below. \\

In each series, the variance is much larger than the mean, which suggests that the data exhibit overdispersion. There is also positive excess skewness, which means that each count is skewed towards the right; this is not surprising since there are a large concentration of zeroes for each type of breach. Furthermore, all the series exhibit excess kurtosis, so they have relatively fatter tails compared to the normal distribution. \\

Univariate Analysis

The INAR model of Section 4.2.4. is the basic dynamic model for a series of count data. However, it implies a marginal Poisson distribution, for which the mean is equal to the variance. Therefore, it is not compatible with data featuring overdispersion, that is, where the variance is much larger than the mean, as seen in Table 2. Consequently, the risk can be underestimated in such a framework. The INAR model can be extended for more flexibility by introducing stochastic intensity. This leads to the univariate Binomial Autoregressive Process (NBAR) [see Gourieroux, Lu (2019)]. We present the univariate model in this section which will be used in this analysis (the bivariate NBAR will be presented in Section 7.3.1.).

Univariate NBAR

The process is defined by its state space representation with a state variable $X_t$ interpretable as stochastic intensity. This representation is as follows:

In total, there are three parameters to be estimated: $\beta$, $\delta$ and $c$. More specifically, $\beta$ characterizes not only the serial dependence of the process, but also the level of conditional overdispersion. The process described above can be represented by the causal chain:

that is a network with one hidden layer and one hidden neuron. Note that the NBAR process is an affine process with conditional Laplace transform [Gourieroux, Lu (2019), Proposition 2]:

where:

and $\rho = \beta c$ and the sequence $(c_h)$ defined by $c_h=c \frac{1-\rho^h}{1-\rho}$ or equivalently, $\beta c_h = \rho \frac{1-\rho^h}{1-\rho}$. This process is strictly stationary if $\rho<1$ and its stationary distribution is obtained for $h$ tending to infinity. The associated Laplace transform is:

and corresponds to the negative binomial distribution $NB(\delta,\rho)$.

Estimation

The NBAR process is a nonnegative Markov process where the distribution is characterized by its conditional Laplace transform at horizon 1. This transition depends on the parameters $\beta$, $c$, $\delta$, by means of $\rho = \beta c$ and $\delta$ only. Therefore, the parameters $\rho = \beta c$ and $\delta$ are identifiable [Gourieroux and Lu (2019)], but not $\beta$ and $c$ separately. Then, we can consider the following two estimation methods: \\

(1) Linear Regression \\

By the law of iterative expectation, it can be shown that:

Therefore, the parameters $\rho$ and $\rho\delta$ (resp. $\delta$) can be estimated using OLS via a regression of $Y_{t}$ on $Y_{t-1}$. However, this method is not asymptotically efficient due to the presence of conditional heteroscedasticity in the NBAR process.\\

(2) Maximum Likelihood Estimation \\

The likelihood function of the NBAR process is given by [Gourieroux, Lu (2019), eq. 4.1]:

where:

and $\Gamma$ denotes the gamma function. Unlike OLS, the MLE estimator is asymptotically efficient. We present below the estimates for $\delta$ and $\rho$ and their associated standard errors under the OLS and MLE methods.

As expected, the estimates of the parameter $\rho$ are positive, to capture the overdispersion, and smaller than one, that is compatible with the stationarity of the count processes. Moreover, the largest persistences are for HACK and ELET. We also observe that the OLS estimates of the persistence (overdispersion) parameter $\rho$ are always smaller than their ML counterparts. This provides some insight on the finite sample bias due to the omission of the conditional heteroscedasticity in the OLS approach.

Univariate FELD

Proof: See Appendix B.2.1. \\

We now present the results of the FELD using the MLE estimated NBAR models for each of the four series for varying risk aversion $u$ and horizon $h$. We focus on the term:

which captures the marginal effect of one additional breach $Y_t$ on the Total FELD.

In each graph, the marginal effect of $Y_t$ is plotted against the forecast horizon $h$. As $h$ falls, the marginal effect of the history $Y_t$ diminishes exponentially. The darker lines represent higher risk aversion parameter $u$. Thus, higher risk aversion means a higher marginal effect of $Y_t$. Intuitively, this implies that if the firm is more risk averse, an additional breach that occurs today will yield higher uncertainty in the forecast (measured by the marginal effect on the Total FELD). Moreover, the marginal effects seem to be higher for HACK and ELET, and lower for DISC and INSD. As seen in Section 3.4, the decompositions in Figure 13 also have interpretation in terms of spot and forward values of derivatives written on the number of cyber events. This is related to the literature on the pricing of cyberinsurance contracts [Fahrenwaldt et al. (2018)].

Bivariate Analysis

The four series correspond to different types of cyber operational risks. Among them two of the risks correspond to breaches with malicious intent, i.e. HACK and INSD. We will focus on the joint analysis of these two risks\footnote{It is possible to study the four risks together, but with more complex models [see online appendix B3 for the extension to the $K$-dimensional NBAR.]}, which has to account for both cross-sectional and serial dependencies between the series. A first insight on the cross-sectional dependence is through the joint stationary distributions. This is illustrated in Figure 14, where on the left hand side we provide a plot of values $(Y_{1,t},Y_{2,t})$ and on the right hand side a plot of the associated Gaussian ranks ($\pi(Y_{1,t})$,$\pi(Y_{2,t})$), i.e. the estimated copula after the Gaussian transform\footnote{The raw data are first ranked by increasing order. Rank$(Y_{j,t})$, the rank of $Y_{j,t}$, is valued in $[1,...,T]$. Then, $\pi(Y_{j,t})=\Phi^{-1}\left[\text{Rank}(Y_{j,t})/T\right]$, where $\Phi$ is the cumulative distribution function of the standard normal distribution.}\footnote{This Gaussian transformed unconditional copula is completed on the count of occurrences. It differs from a copula completed from losses that account for severities [see Eling, Jung (2018)].}. \\

Such plots are usually given for continuous variables or for count variables taking a large number of values. In our framework, the discrete nature of the variable has to be taken into account in the interpretation. For instance, we observe an overweighting of zero counts for the HACK variable, which can be seen in the left plot of Figure 14, where the joint density rests “against" its axis. The right panel in Figure 14 also shows that the two series features some right tail independence.

Bivariate NBAR

Let us denote $Y_{1t},Y_{2t}$ the two count series. We can now introduce three state variables interpretable as specific and common stochastic intensities, denoted $X_{1,t},X_{2,t}$ and $Z_t$, respectively. The nonlinear state space representation becomes:

The dimensionality and the presence of stochastic intensity factors now lead to 9 parameters to be estimated: $\alpha_1$, $\alpha_2$, $\beta_1$, $\beta_2$, $\delta_1$, $\delta_2$, $\sigma_1$, $\sigma_2$, and $\delta$. The process described above corresponds to a more complicated causal scheme with one hidden layer and three hidden neurons.

It is easily checked that the bivariate NBAR model is an affine model. Its conditional Laplace transform at horizon $1$ is given by [see Appendix A.2.7]:

where:

We see from the expression of the conditional Laplace transform that all parameters are identifiable.

Estimation

The increase of the parameter dimension and the introduction of a common stochastic intensity lead to a larger number of identifiable parameters to be estimated, equal to 9. A first estimation to consider is by applying a Vector Autoregressive (VAR) representation based on the linear prediction formula [see Appendix A.2.8]:

Using OLS to estimate the model $Y_t = C + AY_{t-1} +\varepsilon_t$, we obtain the following estimates: \\

The corresponding eigenvalues of $\hat{A} =

$ are $\lambda_1 = 0.639$ and $0.360$, which suggests that the process is stationary in the conditional mean.\\

More generally, we can apply a Method of Moments (MM) procedure based on the unconditional pairwise moment restrictions of the form:

valued for different $u$ and $v$, by considering the orthogonality between prediction errors and past values. A global estimation of the 9 parameters can be based on the unconditional pairwise moments given in (ref) after selecting at least 9 quadruples $(u_1,u_2,v_1,v_2)$ of linearly independent moment restrictions\footnote{We choose quadruplets which correspond to a range of different risk aversion scenarios. For instance, the quadruplet (0.41,0.01,0.41,0.01) corresponds to the scenario where there is high risk aversion on only on the series $Y_{1,t}$. Likewise, the quadruplet (0.41,0.41,0.01,0.01) reflects the case where there is high risk aversion for both series, but only at time $t$.}. We also include 6 additional moment conditions implied by the OLS first order conditions. Applying a generalized method of moment procedure for the 15 moment conditions, we obtain the following estimates displayed in Table 4. \\

These estimates imply $ \hat{\hat{C}} =

\ \hat{\hat{A}} =

$. Hence, our choice of quadruplets not only captures different risk aversion scenarios, but also produce somewhat agreeable estimates with that of OLS. We also note that the estimated values suggest that the process is stationary\footnote{The stationary conditions are given by: $1-\alpha_1-\sigma_1\beta_1>0$, $1-\alpha_2-\sigma_2\beta_2>0$ and $(1-\alpha_1-\sigma_1\beta_1)(1-\alpha_1-\sigma_1\beta_1)>\sigma_1\sigma_2\beta_1\beta_2$ [Gourieroux, Lu (2019), Proposition 3].}.

FELD for Bivariate NBAR

A closed form expression for the conditional Laplace transform at horizon $h$ is difficult to obtain in closed form. However, the dynamic affine property of the bivariate NBAR means that it can be derived numerically by means of recursion. In particular, we have:

where:

The estimated FELD for the model can still be obtained by applying the recursion in (ref) to generate the terms $\Psi(u,h-k|I_{t+k})$ in (ref) and plugging in the estimated values for the 9 parameters from Table 4. We consider two exercises with varying values of $(u_1,u_2)$ [i.e. the risk aversion parameter for HACK and INSD] and of $(Y_{1,t},Y_{2,t})$ [i.e. the last observed historical value for HACK and INSD] to showcase the FELD for the data.

The first exercise is to see how the risk aversion and the historical counts of the two cyberrisks influence the FELD. We compare four scenarios: [1] Low Risk Aversion and Low Historical Counts. [2] High Risk Aversion and Low Historical Counts. [3] Low Risk Aversion and High Historical Counts. [4] High Risk Aversion and High Historical Counts. The results are presented in Figure 15 below. \\

The graph on the left hand side depicts the total FELD up to horizon 10. We observe two important properties of the model. First, the total FELD is increasing in $u$, since the red and blue lines are lower than the green and yellow ones. Second, the total FELD is increasing in $Y$, since the green and blue lines are higher than the red and yellow ones. Indeed, recall that the FELD is a measure of risk. If a firm were to observe high historical counts $Y$ of cyberattacks, then the risk of future attacks should be higher. Likewise, when observing the same level of historical counts of cyberattacks, there is a higher implied risk for a more risk averse firm. \\

On the right hand side graph, we show the decomposition of the total FELD at horizon 10. For ease of comparison between the four scenarios, we express the FELD as a percentage of the total. We have highlighted in red the contribution of risk in updating between horizon 1 and horizon 2 on the FELD at horizon 10. Although there are four different scenarios, we can see that the percentage contribution is roughly 40% in all the cases. This suggests that if we are equally risk averse on the two types of cyber risk, and observe similar levels of historical counts, then our decomposition of the total risk is insensitive to the various levels. \\

Of course, the assumption that the two series have the same risk aversion parameters or the same observed historical values is unrealistic. Hence, we consider a second exercise where these values for HACK and INSD. We compare four new scenarios: [1] Low/High Historical Count for HACK/INSD, High/Low Risk Aversion for HACK/INSD. [2] High/Low Historical Count for HACK/INSD, High/Low Risk Aversion for HACK/INSD. [3] Low/High Historical Count for HACK/INSD, Low/High Risk Aversion for HACK/INSD. [4] High/Low Historical Count for HACK/INSD, Low/High Risk Aversion for HACK/INSD. The results are presented in Figure 16.\\

On the left hand side graph of Figure 16, the yellow and green lines correspond to a high observed historical value for HACK, but a low historical value observed for INSD. These lines are much higher than the blue and red lines, which correspond to the opposite case. Intuitively, this means that observing a high count of breaches due to hacking from an outside party will generate more uncertainty than observing a high count of insider breaches. This is a compelling result since for a firm, it is easier to diagnose the security leak and intitiate preventive measures for future insider breaches than it is for hacking from an outside party. A more straightforward observation is that for high levels of HACK, a lower risk aversion parameter means less uncertainty [i.e. the yellow line is higher than the green one]. \\

On the right hand side graph of Figure 16, we show the decomposition of the total FELD at horizon 10 for this exercise. A notable difference is that now, the decomposition is quite different for the four scenarios. In particular, when updating from horizon 1 to 2 (outlined in red), the scenarios have different contributions to the FELD at horizon 10. When the observed historical count for HACK is low and INSD is high (that is, the first and third bars in the graph), the risk is front loaded since the risk of updating between horizon 1 and 2 accounts for over 40% of the FELD in these cases. On the other hand, when a high historical count of HACK and a low count of INSD is observed (that is, the second and fourth bars in the graph), the risk is has more spread across future horizons. This means that a high count of observed insider breaches implies a more front loaded risk.\\

The exercises above demonstrate two important factors. Firstly, by considering two series together in a bivariate model, we are able to take advantage of the cross-sectional dependencies between the two types of cyber attacks, even if our decomposition measure is a dynamic separation of effects. Secondly, the FELD we have proposed in this paper has as many decompositions as there are $u$ and $Y_t$. In particular, we see that in exercise 1 and exercise 2, the decompositions can be very different and dependent on risk aversion and observed historical counts. For an insurance firm that specializes in cyber risk, these types of scenario analysis can help price their insurance products, or allocate resources efficiently to hedge against future cyber attacks for a range of customers.

Concluding Remarks

In this paper, we have introduced decomposition formulas for the analysis of global forecast errors in nonlinear dynamic models. These formulas are based on functional measures of nonlinear forecast with respect to either their transition densities in the Forecast Error Kullback Decomposition, or the conditional log-Laplace transform in the Forecast Error Laplace Decomposition. These measures and their decompositions are extensions of the Forecast Error Variance Decomposition (FEVD) used in the linear dynamic framework. This latter decomposition is for global shocks on the current value of the process of interest, without trying to define (identify) the sources of the global shocks. The advantage is that this measure and its decomposition with respect to horizon and updating are identifiable, as are also the FEKD and FELD. Such decompositions could be extended to functional measures of economic or financial interest such as conditional Lorenz curves used in inequality analysis or conditional quantiles [Montes-Rojas (2019)]. This is left for future research.

References

Al-Osh, M. and A., Azaid (1987): “First-Order Integer Valued Autoregressive (INAR(1)) Processes", Journal of Time Series Analysis, 8, 261-275. \\

Arrow, K. (1965): “Aspects of the Theory of Risk Bearing", The Theory of Risk Aversion, Helsinki, Reprinted in “Essays in the Theory of Risk Bearing", Markham Publishing Co., Chicago, 1971, 90-109.\\

Barati, M. and B., Yankson (2022): “Predicting the Occurrence of Data Breach," International Journal of Information Management Data Insights, 2. \\

Cebula, J., and L. Young (2010): “A Taxonomy of Operational Cyber Security Risks", Technical Note, Software Engineering Institute, Carnegie Mellon University.\\

Chauvet, M., and S. Potter (2005): “Forecasting Recession Using the Yield Curve", Journal of Forecasting, 24, 77-103. \\

Cox, J., Ingersoll, J., and S., Ross (1985): “A Theory of the Term Structure of Interest Rates", Econometrica, 53, 385-407. \\

Cuchiero, C., Filipovic, D., Mayerhofer, E. and J., Teichmann (2011): “Affine Processes on Positive Semi Definite Matrices", The Annals of Applied Probability, 21, 397-463. \\

Darolles, S., Gourieroux, C., and J., Jasiak. (2006): “Structural Laplace Transform and Compound Autoregressive Models". Journal of Time Series Analysis, 27, 477-503. \\

Doob, J. (1953): “Stochastic Processes", Wiley and Sons. \\

Duffie, D., Filipovic, D., and W., Schachermayer. (2003): “Affine Processes and Applications in Finance". The Annals of Applied Probability, 13, 984-1053. \\

Eling, M., Ibragimov, R., and D., Ning, (2023): “Time Dynamics of Cyber Risk", University of St. Gallen. \\

Eling, M., and K., Jung (2018): “Copula Approaches for Modelling Cross-Sectional Dependence of Data Breach Losses", Insurance: Mathematics and Economics, 82, 167-180. \\

Eling, M., and N. Loperfido (2017): “Data Breaches: Goodness of Fit Pricing and Risk Management", Insurance: Mathematics and Economics, 75, 126 - 136. \\

Embrechts, P. (2000): “Actuarial versus Financial Pricing of Insurance", J. Risk Finance, 1, 17-26. \\

Estrella, A., and F. Mishkin (1998): “Predicting US Recessions: Financial Variables as Leading Indicators", Review of Economics and Statistics, 80, 45-61. \\

Fahrenwaldt, M., Weber, S., and K., Weske (2018): “Pricing of Cyber Insurance Contracts in a Network Model", ASTIN Bulletin, 48, 1175-1218. \\

Feller, W. (1971). “An Introduction to Probability Theory and Its Applications", (2nd Edition), Vol 2, Wiley, New York. \\

Fishburn, R., and R., Vickson (1978): “Theoretical Foundations of Stochastic Dominance", in Stochastic Dominance: An Approach to Decision Making Under Risk, Whitmore, G. and Findlay, M. (eds). DC Health, London.\\

Gourieroux, C., and J., Jasiak (2006): “Autoregressive Gamma Processes", Journal of Forecasting, 25, 129-152. \\

Gourieroux, C., Jasiak, J., and R., Sufana (2009): “The Wishart Autoregressive Process of Multivariate Stochastic Volatility", Journal of Econometrics, 150, 167-181. \\

Isakin, M., and P., Ngo (2020): “Variance Decomposition Analysis for Nonlinear Economic Models". Oxford Bulletin of Economics and Statistics, 82, 1362-1374.\\

Kauppi, H., and P., Saikkonen (2008): “Predicting US Recessions with Dynamic Binary Response Models", Review of Economics and Statistics, 90, 777-796. \\

Lanne, M., and H., Nyberg. (2016): “Generalized Forecast Error Variance Decomposition for Linear and Nonlinear Multivariate Models". Oxford Bulletin of Economics and Statistics, 78(4), 595-603.\\

Lu, Y., Zhang, J., and W., Zhu (2024): “Cyber Risk Modelling: A Discrete Multivariate Count Process Approach," forthcoming, Scandinavian Actuarial Journal.\\

Markowitz, H. (1952): “Portfolio Selection", The Journal of Finance, 7,77-91. \\

Markowitz, H. (2000): “Mean-Variance Analysis in Porfolio Choice and Capital Markets", Wiley and Sons, Vol. 66. \\

McKenzie, E. (1985): “Some Simple Models for Discrete Variate Time Series", Water Resources Bulletin, 21, 645-650. \\

Montes-Rojas, G. (2019): “Multivariate Quantile Impulse Response Functions", Journal of Time Series Analysis, 40, 730-752. \\

Muirhead, R. (1982): “Aspects of Multivariate Statistical Theory", Wiley and Sons. \\

Pratt, J. (1964): “Risk Aversion in the Small and in the Large", Econometrica, 32, 129-136. \\

Rothschild, M., and J., Stiglitz (1970): “Increasing Risk: I. A Definition", Journal of Economic Theory, 2, 225-243. \\

Schwartz, G. and S. Sastry (2014): “Cyber-insurance framework for large scale interdependent networks", in Proceedings of the 3rd International Conference on High Confidence Networked Systems, 145-154, New-York, The Association on Computing Machinery. \\

Sun, H., Xu. M., and P., Zhao (2021): “Modelling Malicious Hacking Data Breach Risks", North American Actuarial Journal, 25, 484-502. \\

Vickson, R. (1975): “Stochastic Dominance for Decreasing Absolute Risk Aversion", Journal of Financial and Quantitative Analysis, 10, 799-811. \\

Whittle, P. (1963): “Prediction and Regulation", Bowman, new edition in 1983. \\