# Convergence of Computed Dynamic Models with Unbounded Shock

Kenichiro<EMAIL_ADDRESS>[ Kōsaku
<EMAIL_ADDRESS>[ Department of Statistical
Science, Fox School of Business, Temple University  Center for Advanced
Intelligence Project, Riken

###### Abstract

This paper studies the asymptotic convergence of computed dynamic models when
the shock is unbounded. Most dynamic economic models lack a closed-form
solution. As such, approximate solutions by numerical methods are utilized.
Since the researcher cannot directly evaluate the exact policy function and
the associated exact likelihood, it is imperative that the approximate
likelihood asymptotically converges– as well as to know the conditions of
convergence– to the exact likelihood, in order to justify and validate its
usage. In this regard, [Fernandez-Villaverde, Rubio-Ramirez, and Santos
(2006)] show convergence of the likelihood, when the shock has compact
support. However, compact support implies that the shock is bounded, which is
not an assumption met in most dynamic economic models, e.g., with normally
distributed shocks. This paper provides theoretical justification for most
dynamic models used in the literature by showing the conditions for
convergence of the approximate invariant measure obtained from numerical
simulations to the exact invariant measure, thus providing the conditions for
convergence of the likelihood.

Dynamic economic models,

Convergence,

Computation,

###### keywords:

##  1 Introduction

This paper studies the convergence of dynamic economic models. While dynamic
economic models have become a central tool for research and policy, most do
not have a closed-form solution. Due to this, the policy function of these
economic models are approximated by numerical methods. This approximation
means that the researcher can only evaluate the approximated transition
function associated with the approximated invariant measure, rather than the
exact invariant measure implied by the exact transition function. Given that
the researcher cannot evaluate the exact measure, it is natural to ask whether
the approximate measure converges to the exact measure, at least
asymptotically. If, for example, it does not converge, or the conditions for
convergence are not met in practice, then the validity of the estimated
economic model and its output comes into question. It is, therefore, critical
that there is a theoretical foundation that provides the conditions for
convergence to justify the usage of these dynamic economic models.

As a response, much econometric analysis has been done to provide this
theoretical foundation. For example, [Santos (2004)] and [Santos and Peralta-
Alva (2005)] provide the foundations of simulations of approximate solutions
for stochastic dynamic models by studying its accuracy properties, showing
that the computed moments from the numerically approximated policy converge to
the exact moments as the approximation errors of the computed solutions go to
zero. Further, and more relevant to this paper, [Fernandez-Villaverde, Rubio-
Ramirez, and Santos (2006)] extends the results of [Santos and Peralta-Alva
(2005)] to the convergence of the likelihood of computed economic models,
providing conditions for which the approximated likelihood functions converges
to the exact likelihood. While the convergence results in [Fernandez-
Villaverde, Rubio-Ramirez, and Santos (2006)] provide some justification for
dynamic economic models, one assumption it employs to obtain their result is
rarely met. This assumption is the compactness of the state space, which
implies that the support of the shock of a dynamical system is bounded.
Although this assumption is standard in the numerical literature, it excludes–
among others– dynamical models with normally distributed shocks. As assuming a
normally distributed shock is standard in empirical studies, in which the
evaluation of the likelihood is done by the usage of the Kalman filter (Smets
and Wouters, 2007), it is simply vital that the results in Fernandez-
Villaverde, Rubio-Ramirez, and Santos (2006) extend to non-compact support,
i.e., unbounded shock. For example, recent works by Stachurski (2002),
Nishimura and Stachurski (2005), Kamihigashi (2007) and Kamihigashi and
Stachurski (2016) study the asymptotic invariant measure of the stochastic
neoclassical growth model without compactness of the shocks and states. The
purpose of this paper is to relax the compactness assumption for the
convergence of the approximated invariant measure, providing the theoretical
foundation and justification for these models.

The rest of this paper is organized as follows. Section 2 gives the set-up of
dynamic economic models and preliminary of the Markov operator. Section 3
presents our result on the convergence of the invariant measure. In Section 4,
we derive error bounds for these approximations. Section 5 presents our main
result on the convergence of computed likelihoods.

##  2 Model Set-Up and Preliminaries

We follow the set of notations and models in Santos and Peralta-Alva (2005).
The equilibrium law of motion of the state variables is specified by a
dynamical system of the form

$\displaystyle s_{n+1}$
$\displaystyle=\varphi\left(s_{n},\varepsilon_{n+1}\right),\qquad
n=0,1,2,\ldots,$

where $s_{n}$ is a vector of state variables that characterize the evolution
of the system. The vector $s_{n}$ belongs to a measurable state space
$\left(S,\mathcal{S}\right)$. We endow $S$ with its relative Borel
$\sigma$-algebra, which we denote by $\mathcal{S}$. The variable $\varepsilon$
is an independent and identically distributed shock, which is defined on the
sample space $\left(E,\mathcal{E}\right)$. The distribution of the shock
$\varepsilon$ is given by a stochastic kernel $Q:S\times
E\rightarrow\left[0,1\right]$, where $Q\left(s,A\right)$ is the probability of
realizing the event $A\in\mathcal{E}$, given that the current state is $s\in
S$.

Given a random dynamical system, one can define a transition probability on
the state space in the following way. Define the transition probability
function as

$\displaystyle P\left(s,A\right)$
$\displaystyle=Q\left(\left\\{\varepsilon\mid\varphi\left(s,\varepsilon\right)\in
A\right\\}\right).$

The transition function, $P:S\times\mathcal{S}\rightarrow\left[0,1\right]$, is
defined by

$\displaystyle P\left(s,A\right)$
$\displaystyle=Q\left(s,\varphi^{-1}\left(A\right)\right).$

Let $B\left(S\right)$ be the set of all bounded $S$-measurable real valued
functions on $S$, with sup norm
$\left|f\right|=\sup_{S}\left|f\left(s\right)\right|$. The Markov operator
associated with $P$ is defined as

$\displaystyle Tf\left(s\right)$ $\displaystyle\triangleq\int
f\left(t\right)P\left(s,dt\right)$ (1) $\displaystyle=\int
f\left(\varphi\left(s,\varepsilon\right)\right)Q\left(s,d\varepsilon\right)$

For any given initial condition $\mu_{0}$ on $\mathcal{S}$, the evolution of
future probabilities, $\left\\{\mu_{n}\right\\}$, can be specified by the
following operator $T^{*}$ that takes the space

$\displaystyle\mu_{n+1}$
$\displaystyle=\left(T^{*}\mu_{n}\right)\left(A\right)={\displaystyle\int}P\left(s,A\right)\mu_{n}\left(ds\right),$

for all $A$ in $\mathcal{S}$ and $n\geq 0$. The adjoint $T^{*}$of $T$ is
defined by the formula

$\displaystyle\left(T^{*}\mu\right)\left(A\right)$ $\displaystyle=\int
P\left(t,A\right)\mu\left(dt\right).$

We maintain the following basic assumptions.

######  Assumption 2.1.

The space of sets $S$ and $E$ are both locally compact and $\sigma$-compact.

Locally compact means that for each point $x\in S$, there is some compact
subspace $C$ of $S$ that contains a neighborhood of $x\in S$. Further,
$\sigma$-compact is a countable union of compact spaces. Note that the space
$\mathbb{R}^{d}$ is both locally compact and $\sigma$-compact. A space that is
both locally compact and $\sigma$-compact can be written as an increasing
union of countably many open sets, each of which is compact and closed. In
Santos and Peralta-Alva (2005), they impose the compactness assumption on both
states, $S$ and $E$, which, again, is not an assumption met in most dynamic
economic models used in empirical studies. In an important distinction, we
relax this restriction to the non-compact case, which allows us to use the
whole Euclidean state, $S=\mathbb{R}^{d}$, and unbounded distributions, such
as the normal distribution.

Recall that the probability measure $P$ is called tight if for all
$\epsilon>0$ there is a compact set $K\subset\mathcal{S}$ such that
$P\left(K\right)\geqq 1-\epsilon$. Any probability measure on the complete
separable metric space is tight.

######  Assumption 2.2.

The Markov operator $T^{*}$ has a unique fixed point $\mu_{0}$:
$T^{*}\mu_{0}=\mu_{0}$.

A sufficient condition for Assumption 2.2 is that there exists a point,
$s_{0}\in S$, such that, for any point $s\in S$, any neighborhood $U$ of
$s_{0}$ and any integer $k\geq 1$, we have $P^{nk}\left(s,U\right)>0$ (see
Futia, 1982, Section 3.2).

######  Assumption 2.3.

Function $\varphi{:}\ S\times E\rightarrow S$ is jointly measurable. Moreover,
for every continuous function $f{:}\ S\rightarrow\mathbb{R}$,

$\displaystyle{\displaystyle\int}f\left(\varphi\left(s_{j},\varepsilon\right)\right)Q\left(d\varepsilon\right)$
$\displaystyle\underset{}{\rightarrow}{\displaystyle\int}f\left(\varphi\left(s,\varepsilon\right)\right)Q\left(d\varepsilon\right)\;as\;s_{j}\rightarrow
s.$

Assumption 2.3 is the same as Assumption 2 in Santos and Peralta-Alva (2005).

In most cases, the researcher does not know the exact form of the transition
equation $\varphi$, and only has access to the numerical approximation of the
transition equation, $\varphi_{j}$, with index $j$. The index $j$ indicates
the approximation and implies that, as $j$ goes to infinity, the
approximation, $\varphi_{j}$, converges to the exact value (the metric of
convergence is defined later). Every numerical approximation $\varphi_{j}$
defines the transition probability $P_{j}$ on $\left(S,\mathcal{S}\right)$.
Given an approximation $\varphi_{j}$, we define the corresponding
approximation of the transition probability as

$\displaystyle P_{j}\left(s,A\right)$
$\displaystyle=Q\left(\left\\{\varepsilon\mid\varphi_{j}\left(s,\varepsilon\right)\in
A\right\\}\right),$

and define the approximated transition function,
$P_{j}:S\times\mathcal{S}\rightarrow\left[0,1\right]$, as

$\displaystyle P_{j}\left(s,A\right)$
$\displaystyle=Q\left(s,\varphi_{j}^{-1}\left(A\right)\right).$

The Markov operator associated with $P_{j}$ is defined as

$\displaystyle T_{j}f\left(s\right)$ $\displaystyle\triangleq\int
f\left(t\right)P_{j}\left(s,dt\right)$ (2) $\displaystyle=\int
f\left(\varphi_{j}\left(s,\varepsilon\right)\right)Q\left(s,d\varepsilon\right).$

The evolution of future probabilities, $\left\\{\mu_{n}^{j}\right\\}$, can be
specified by the following operator $T_{j}^{*}$ that takes the space

$\displaystyle\mu_{n+1}^{j}$
$\displaystyle=\left(T_{j}^{*}\mu_{n}^{j}\right)\left(A\right)={\displaystyle\int}P_{j}\left(s,A\right)\mu_{n}^{j}\left(ds\right),$

for all $A$ in $\mathcal{S}$ and $n\geq 0$. The adjoint $T_{j}^{*}$ of $T_{j}$
is defined by

$\displaystyle\left(T_{j}^{*}\mu\right)\left(A\right)$ $\displaystyle=\int
P_{j}\left(t,A\right)\mu\left(dt\right).$

Every numerical approximation $\varphi_{j}$ satisfies a structure parallel to
that of the above Assumptions 2.2 and 2.3. We further assume:

######  Assumption 2.4.

The Markov operator $T_{j}^{*}$ has a unique fixed point $\mu_{0}^{j}$:
$T_{j}^{*}\mu_{0}^{j}=\mu_{0}^{j}$, for all $j$.

######  Assumption 2.5.

For each $j$, the function $\varphi_{j}{:}\ S\times E\rightarrow S$ is jointly
measurable. Moreover, for every continuous function $f{:}\
S\rightarrow\mathbb{R}$,

$\displaystyle{\displaystyle\int}f\left(\varphi_{j}\left(s_{j},\varepsilon\right)\right)Q\left(d\varepsilon\right)$
$\displaystyle\underset{}{\rightarrow}{\displaystyle\int}f\left(\varphi_{j}\left(s,\varepsilon\right)\right)Q\left(d\varepsilon\right)\;as\;s_{j}\rightarrow
s.$

##  3 Convergence of the Invariant Distribution

Now, recall the convergence of probability measures on $S$. When the state
space $S$ is separable, we can introduce a metric $D$ in the space of
probability measures on $S$, such that $\lim_{n}D\left(\mu_{n},\mu\right)=0$
if and only if $\mu_{n}$ converges in law to $\mu$. Specifically, the metric
we use is the Fortet-Mourier metric (Dudley, 2002, Section11.3):

$\displaystyle D\left(\mu_{n},\mu\right)$ $\displaystyle=\sup_{f\in
BL\left(S\right)}\left|\int_{S}f\left(s\right)d\mu_{n}-\int_{S}f\left(s\right)d\mu\right|,$
(3)

where the supremum $\sup_{f\in BL\left(S\right)}$ is taken over all bounded
Lipschitz continuous functions defined on $S$: $BL\left(S\right)$.

The main question we answer in this paper is the following: How strong of a
topology is sufficient for the approximate transition equation, $\varphi_{j}$,
to converge to the true transition equation, $\varphi$, in order for the
approximate invariant measure, $\mu_{n}$, to converge to the convergence in
law distance eq. (3). Santos and Peralta-Alva (2005), assuming that the state-
space is compact, showed that convergence under the following topology is
sufficient to prove the convergence of the invariant measure. Endow the metric
in the space of functions $\varphi$ and $\hat{\varphi}$ as

$\displaystyle\max_{s\in
S}\left[{\displaystyle\int}\left\|\varphi\left(s,\varepsilon\right)-\hat{\varphi}\left(s,\varepsilon\right)\right\|Q\left(d\varepsilon\right)\right]$

where $\left\|\cdot\right\|$ is the max norm in $\mathbb{R}^{l}$. This metric
only works under the compactness assumption on $S$. To consider the functional
approximation of the transition equation, $\varphi$, under non-compactness,
this uniform topology is not practical. In the following, we extend the state-
space, $S$, to non-compactness and weaken the uniform convergence topology of
the functional approximation to a local uniform topology.

First, note that $BL\left(S\right)$ can be relaxed to infinitely continuously
differentiable functions on $S$: $C^{\infty}\left(S\right)$ by using the
mollifier method. Then, we have the following lemma:

######  Lemma 3.1.

The limit, $\lim_{n}D\left(\mu,\mu_{n}\right)=0$, holds if and only if

$\displaystyle\lim_{n\rightarrow\infty}\sup_{f\in
C^{\infty}\left(S\right)}\left|\int_{S}f\left(s\right)d\mu_{n}-\int_{S}f\left(s\right)d\mu\right|$
$\displaystyle=0.$

Note that by Assumptions 2.2 and 2.3, each $\varphi_{n}$ defines the
associated pair $\left(P_{j},T_{j}\right)$. The adjoint $T_{j}^{*}$ of $T_{j}$
is

$\displaystyle\left\langle T_{j}f,\mu_{j}\right\rangle$
$\displaystyle=\int\int
f\left(\varphi_{j}\left(s,\varepsilon\right)\right)Q\left(s,d\varepsilon\right)d\mu_{j}\left(s\right)$
$\displaystyle\left\langle f,T_{j}^{*}\mu_{j}\right\rangle$
$\displaystyle=\int\int
f\left(\varphi_{j}\left(s,\varepsilon\right)\right)Q\left(s,d\varepsilon\right)d\mu_{j}\left(s\right).$

Moreover, there always exists an invariant distribution
$\mu_{j}^{*}=T_{j}^{*}\mu_{j}^{*}$.

Given Lemma 3.1, we have the following result:

######  Proposition 3.1.

Suppose Assumptions 2.1, 2.2, and 2.3 are satisfied for each approximated
model $\varphi_{i}$ and $\varphi_{0}$. Then, a sufficient condition for a
sequence of the measure $\mu_{j}$, associated with $T_{j}$, to converge to
$\mu_{0}$, associated with $T_{0}$, in the sense of eq. (3), is the strong
convergence of $T_{j}$ to $T_{0}$.

Since $S$ is completely regular, it has Stone-Cech compactification
$\beta\left(S\right)$:

######  Theorem 3.1.

(Stone-Cech compactification: Munkres, 2000, Theorem 38.2). Let $S$ be a
completely regular space. Then, there exists a compactification
$\beta\left(S\right)$ of $S$ having the property that every bounded continuous
function $f{:}\ S\rightarrow\mathbb{R}$ extends uniquely to a continuous
function of $\beta\left(S\right)$ into $\mathbb{R}$.

The Stone-Cech compactification, $\beta\left(S\right)$, of which $S$ is a
dense subspace, satisfies the property that each bounded continuous function,
$f{:}\ S\rightarrow\mathbb{R}$, has a continuous extension,
$g:\beta\left(S\right)\rightarrow\mathbb{R}$. We endow the metric in the space
of functions defined on the locally compact and $\sigma$-compact space, $S$.
For any two vector-value functions $\varphi$ and $\hat{\varphi}$, let
$d\left(\cdot,\cdot\right)$ be

$\displaystyle d\left(\varphi,\hat{\varphi}\right)$ $\displaystyle=\max_{f\in
C^{\infty}\left(\beta\left(S\right)\right)}\max_{s\in\beta\left(S\right)}\left[{\displaystyle\int}\left|f\left(\varphi\left(s,\varepsilon\right)\right)-f\left(\hat{\varphi}\left(s,\varepsilon\right)\right)\right|Q\left(d\varepsilon\right)\right].$
(4)

The metric in eq. (4) is weaker than the metric of Santos and Peralta-Alva
(2005), and extends to the non-compact state space. In this section,
convergence of the sequence of functions $\left\\{\varphi_{j}\right\\}$ is in
this distance, as this metric can accommodate the noncontinuous functions
$\varphi$ and $\hat{\varphi}$. Although we will impose continuous
differentiability on $\varphi$ for the convergence of the approximate
likelihood studied in Section 5, the metric $d\left(\cdot,\cdot\right)$ is
sufficient to guarantee the convergence of the invariant distribution.

Then, we have the following theorem:

######  Theorem 3.2.

Let $\left\\{\varphi_{j}\right\\}$ be a sequence of functions that converge to
$\varphi$, in the sense of $d\left(\cdot,\cdot\right)$ in eq. (4). Let
$\left\\{\mu_{j}^{*}\right\\}$ be a sequence of probabilities on $\mathcal{S}$
associated with $\left\\{\varphi_{j}\right\\}$, such that
$\mu_{j}^{*}=T_{j}^{*}\mu_{j}^{*}$, for each $j$. Under Assumptions 2.1 and
2.2, if $\mu^{*}$ is a weak limit point of $\left\\{\mu_{j}^{*}\right\\}$,
then $\mu^{*}=T^{*}\mu^{*}$.

This theorem asserts the bilinear convergence of $T_{j}^{*}\mu_{j}^{*}$ to
$T^{*}\mu^{*}$ in the weak topology.

##  4 Error Bounds

In this section, we study the error bounds of these approximations under non-
compactness. The error bounds are important for two reasons. First, in
numerical applications, it is often desirable to bound the size of the
approximation error in order to know the theoretical limit of the
approximation. Second, computations cannot go on forever and must stop in
finite time. Hence, knowing the error bounds can dictate an efficient stopping
criteria to minimize computational cost while ensuring convergence. As such,
Santos and Peralta-Alva (2005) give a bound on the size of the approximation
error under the compactness assumption.

To begin, we introduce the notion of compactness for the Markov operator. The
Markov operator $T$ is compact if the image $T\left(bX\right)$ has compact
closure in $X$, where $bX=\left\\{x\in X|\left\|x\right\|\leq 1\right\\}$. The
Markov operator $T$ is quasi-compact if there is a unique compact operator $L$
and an integer $n$ such that

$\sup_{x\in bX}\left\|T^{n}x-Lx\right\|<1.$

If the above quasi-compactness is satisfied, one can obtain the convergence of
the sequence of operators, $\left\\{T^{n}\right\\}$, to the invariant
probability at a geometric rate. The following theorem gives this result.

######  Theorem 4.1.

(Yosida and Kakutani, 1941). Let $T$ be a stable, quasi-compact Markov
operator defined by eq. (1) satisfying Assumptions 2.1, 2.2, and 2.3. Then,
there exist constants, $C$, with $\varepsilon>0$, such that

$\displaystyle\sup_{s\in\beta\left(S\right)}\left\|T^{n}f\left(s\right)-T^{*}f\left(s\right)\right\|$
$\displaystyle\leq\frac{C}{\left(1+\varepsilon\right)^{n}}.$

The following theorem bounds the approximation error between the expected
values of $f$ over the true invariant measure $\mu^{*}$ and the approximate
invariant measure $\hat{\mu}^{*}$ of $\hat{\varphi}$.

######  Proposition 4.1.

Let $f$ be a Lipschitz function with constant $L$. Suppose we have a numerical
approximation $\hat{\varphi}$ with
$d\left(\hat{\varphi},\varphi\right)\leq\delta$, for some $\delta>0$. Then,

$\left|\int f\left(s\right)\mu^{*}\left(ds\right)-\int
f\left(s\right)\hat{\mu}^{*}\left(ds\right)\right|\leq\frac{Ld\left(\hat{\varphi},\varphi\right)}{\varepsilon},$

where $\mu^{*}$ is the unique invariant measure of the exact $\varphi$, and
$\hat{\mu}^{*}$ is a unique invariant measure of $\hat{\varphi}$.

Quasi-compact operators enjoy a very useful property in Theorem 4.1.
Furthermore, quasi-compact operators are easily recognizable. In fact, we find
that most operators are quasi-compact (see Futia, 1982). In our circumstance,
Assumption 2.2 guarantees the quasi-compactness of the Markov operator.

##  5 Application: Convergence of Computed Likelihood

Given the convergence of the invariant measure in the previous section, we
prove the convergence of the approximate likelihood in Fernandez-Villaverde,
Rubio-Ramirez, and Santos (2006). We relax the compactness assumption in the
state-space and shock, and prove equivalent results to Fernandez-Villaverde,
Rubio-Ramirez, and Santos (2006), justifying the construction of the
likelihood via the Kalman filter, among others.

###  5.1 Likelihood Induced by Random Dynamical Systems

The equilibrium law of motion of the state space system can be specified as

$\displaystyle s_{t}$ $\displaystyle=\varphi(s_{t-1},\varepsilon_{t};\theta),$
(5) $\displaystyle y_{t}$ $\displaystyle=g(s_{t},\eta_{t};\theta),$ (6)

where eq. (5) is the transition equation, and eq. (6) is the measurement
equation. Here, the variables $\varepsilon_{t}$ and $\eta_{t}$ are tight
random elements and are independent and identically distributed shocks with
values in some Euclidean space, with bounded and continuous densities. Their
distribution is given by the probability measure, $Q$, defined on a measurable
space,
$\left(E,\mathcal{E}\right)\subset\left(\mathbb{R}^{d},\mathcal{B}\left(\mathbb{R}^{d}\right)\right)$.
We do not impose the compactness on the support of $Q$ but impose tightness,
in order to deal with unbounded shocks, such as normally distributed shocks.
The parameter, $\theta\in\Theta\subset\mathbb{R}^{n}$, is a vector of
structural parameters and $\Theta$ is on a compact set. The vector, $y_{t}$,
is the observables in each period, $t$. Let
$\mathcal{Y}_{T}=\left\\{y_{t}\right\\}_{t=1}^{T}$ with
$Y^{0}=\left\\{\emptyset\right\\}$. To avoid singularity, we impose
$\textrm{dim}(\varepsilon_{t})+\textrm{dim}(\eta_{t})\geq\textrm{dim}(Y_{t})$.
And we partition $\left\\{\varepsilon_{t}\right\\}$ into
$\varepsilon_{t}=(\varepsilon_{1,t},\varepsilon_{2,t})$, such that
$\textrm{dim}(\varepsilon_{2,t})+\textrm{dim}(\eta_{t})=\textrm{dim}(y_{t}).$
As in the previous section, we index the approximations by $j$, the numerical
approximation to the transition equations is $\varphi_{j}$, and the
measurement equations is $g_{j}$.

As with the previous section, we assume that each state-space system has an
invariant measure and that invariance measure is absolute continuous with
regard to a Lebesgue measure:

######  Assumption 5.1.

For all $\theta$ and all $j$, there exists a unique invariant distribution for
$\mathcal{S}$, $\mu(S;\theta)$, and $\mu_{j}(S;\theta)$, that has a Radon-
Nikodym derivative with respect to the Lebesgue measure.

The exact likelihood is constructed, using the change of variables formula, as
follows. First we assume that the system can solve the error term, exactly.

######  Assumption 5.2.

For all $\theta$ and $t$, the system of equations eq. (5) and eq. (6) have a
unique solution,

$\displaystyle\eta_{t}$
$\displaystyle=\eta^{t}(\varepsilon_{1}^{t},s_{0},y^{t};\theta),$
$\displaystyle s_{t}$
$\displaystyle=\mathit{s}^{t}(\varepsilon_{1}^{t},s_{0},y^{t};\theta),$
$\displaystyle\varepsilon_{2,t}$
$\displaystyle=\varepsilon_{2}^{t}(\varepsilon_{1}^{t},s_{0},y^{t};\theta),$

and we can evaluate $p(\mathit{v}^{t}(W_{1}^{t},S_{0},y^{t};\theta);\theta)$
and $p(\mathit{w_{2}}^{t}(W_{1}^{t},S_{0},y^{t};\theta);\theta)$ for all
$S_{0}$, $W_{1}^{t},$ and $t$.

We further assume that the observation equation (6) is continuously
differentiable.

######  Assumption 5.3.

For all $\theta$, function $g(\cdot,\cdot,\theta)$ is continuously
differentiable, with bounded partial derivatives.

From Assumptions 5.3, 5.1, and 5.2, we construct the likelihood function by
the change of variable formula:

$p(y_{t}\mid
W_{1}^{t},S_{0},y^{t-1};\theta)=p(\mathit{v}_{t};\theta)p(\mathit{w}_{2,t};\theta)\mid
dy(\mathit{v}_{t},\mathit{w}_{2,t};\theta)\mid,$ (7)

where

$\mid
dy(\mathit{v}_{t},\mathit{w}_{2,t};\theta)\mid=\textrm{det}\left[\begin{array}[]{cc}\nabla\frac{\partial
g}{\partial v_{t}}&\nabla\frac{\partial g}{\partial
w_{2,t}}\end{array}\right].$

Further, we have the following assumption (Fernandez-Villaverde, Rubio-
Ramirez, and Santos, 2006, Assumption 4).

######  Assumption 5.4.

For all $\theta$ and $t$, the model gives some positive probability to the
data $y^{T}$, that is, $p(y_{t}\mid W_{1}^{t},S_{0},y^{t-1};\theta)>\xi\geq 0$
for all $S_{0}$ and $W_{1}^{t}.$

From Assumption 5.4, the likelihood is as follows,

$\displaystyle L(y^{T};\gamma)$
$\displaystyle={\displaystyle\prod_{t=1}^{T}p(y_{t}\mid y^{t-1};\theta)}$
$\displaystyle={\displaystyle\prod_{t=1}^{T}\int\int p(y_{t}\mid
W_{1}^{t},S_{0},y^{t-1};\theta)p(W_{1}^{t},S_{0}\mid
y^{t-1};\theta)dW_{1}^{t}dS_{0}}$
$\displaystyle=\int\left(\int\prod_{t=1}^{T}p(W_{1}^{t};\theta)p(y_{t}\mid
W_{1}^{t},S_{0},y^{t-1};\theta)dW_{1}^{t}\right)\mu^{*}(dS_{0};\theta).$

Next, we also assume that, also for the approximate state-space functions,
$\left\\{\varphi_{j}\right\\}$, and measurement function,
$\left\\{g_{j}\right\\}$, the system can solve the error term, exactly.

######  Assumption 5.5.

For all $j$, the system of equations

$\displaystyle S_{1}$
$\displaystyle=\varphi_{j}(S_{0},(W_{1,1},W_{2,1});\theta),$ $\displaystyle
y_{m}$ $\displaystyle=g_{j}(S_{m},V_{m};\theta)\quad\textrm{for}\quad
m=1,2,\dots,t,$ $\displaystyle S_{m}$
$\displaystyle=\varphi_{j}(S_{m-1},(W_{1,m},W_{2,m});\theta)\quad\textrm{for}\quad
m=2,3,\dots,t,$

has a unique solution,

$\displaystyle V_{j,t}$
$\displaystyle=\mathit{v}_{j}^{t}(W_{1}^{t},S_{0},y^{t};\theta),$
$\displaystyle S_{j,t}$
$\displaystyle=\mathit{s}_{j}^{t}(W_{1}^{t},S_{0},y^{t};\theta),$
$\displaystyle W_{j,2,t}$
$\displaystyle=\mathit{w}_{j,2}^{t}(W_{1}^{t},S_{0},y^{t};\theta),$

and we can evaluate
$p(\mathit{v}_{j}^{t}(W_{1}^{t},S_{0},y^{t};\theta);\theta)$ and
$p(\mathit{w_{j,2}}^{t}(W_{1}^{t},S_{0},y^{t};\theta);\theta)$ for all
$S_{0}$, $W_{1}^{t},$ and $t$.

We also assume that the measurement function, $\left\\{g_{j}\right\\}$, is
continuously differentiable.

######  Assumption 5.6.

For all $j$, functions $g_{j}(\cdot,\cdot,\theta)$ are continuously
differentiable at all points except at a finite number of points. At the
points of differentiability, all partial derivatives are bounded, and the
bounds are independent of $j$.

Then, $dy_{j}(v_{j,t},w_{j,2};\theta)$ exists for all but a finite set of
$S_{0}$ and $W_{1}^{t}$, we have, for all $j$, $\theta$, and $t$,

$p_{j}(y_{t}\mid
W_{1}^{t},S_{0},y^{t-1};\theta)=p(\mathit{v}_{j,t};\theta)p(\mathit{w}_{j,2,t};\theta)\mid
dy_{j}(\mathit{v}_{j,t},\mathit{w}_{j,2,t};\theta)\mid,$

where

$\mid
dy_{j}(\mathit{v}_{j,t},\mathit{w}_{j,2,t};\theta)\mid=\textrm{det}\left[\begin{array}[]{cc}\nabla\frac{\partial
g_{j}}{\partial v_{j,t}}&\nabla\frac{\partial g_{j}}{\partial
w_{j,2,t}}\end{array}\right],$

for all $S_{0}$ and $W_{1}^{t}$, but a finite number of points.

As $j$ goes to infinity, $\varphi_{j}$ and $g_{j}$ converge to their exact
values. Unlike the previous section, convergence of the sequence of functions
$\left\\{\varphi_{j}\right\\}$ and $\left\\{g_{j}\right\\}$ need to be a
stronger topology, which is defined in the following way. For any two vector-
valued functions $\varphi$ and $\hat{\varphi}$, let

$\displaystyle d_{C^{1}}\left(\varphi,\hat{\varphi}\right)$
$\displaystyle=\max_{i\in I}\sup_{s\in
S_{i}}\left[\int\left\|\varphi\left(s,\varepsilon\right)-\hat{\varphi}\left(s,\varepsilon\right)\right\|Q\left(d\varepsilon\right)+\int\left\|\nabla\varphi\left(s\right)-\nabla\hat{\varphi}\left(s\right)\right\|Q\left(d\varepsilon\right)\right],$

where $\left\\{S_{i}\right\\}_{i\in I}$ is an exhaustive sequence of compact
sets of $S$. In this section, convergence of a sequence of functions,
$\left\\{\varphi_{j},g_{j}\right\\}$, should be understood in this norm.

Assumption 5.5 is required for the change of variable formula in eq. (7).
Convergence in $C^{1}$ implies the convergence of the solutions, $v_{j}^{t}$,
$s_{j}^{t},$ and $w_{j,2}^{t}$, from the same argument in the proof of
consistency of Z-estimators. More importantly, the convergence of the
Jacobian,
$\left|dy_{j}\left(\mathit{v}_{j,t},\mathit{w}_{j,2,t};\theta\right)\right|$,
can be derived from the Ascoli-Arzela theorem.

######  Theorem 5.1.

(Ascoli-Arzela). Let $S$ be a separable $\sigma$-compact metric space and
$S_{i}$, $i\in I$, be an exhausting sequence of compact sets. Further, let
$C^{1}(S)$ be the Banach space of complex-valued continuous functions, $f(s)$,
normed by

$\displaystyle\parallel f\parallel$ $\displaystyle=\max_{i\in
I}\left[\sup_{s\in S_{i}}\mid f(s)\mid+\sup_{s\in
S_{i}}\left|f^{\prime}\left(s\right)\right|\right].$

If the following two conditions are satisfied:

$\begin{array}[]{cc}{\displaystyle\sup_{n\geq 1}\sup_{s\in S}\mid
f_{n}(s)\mid<\infty},&\textrm{(equi-bounded)}\\\
{\displaystyle\lim_{\delta\downarrow 0}\sup_{n\geq
1,dis(s^{\prime},s^{\prime\prime})\leq\delta}\mid
f_{n}(s^{\prime})-f_{n}(s^{\prime\prime})\mid=0},&\textrm{(equi-
continuous)}\end{array}$

then the sequence, $\left\\{f_{n}(s)\right\\}\subseteq C(S)$, is relatively
compact in $C(S)$.

Given this, we have the following proposition that proves the convergence of
the approximate likelihood to the exact likelihood with unbounded shock.

######  Proposition 5.1.

Suppose Assumptions 5.3-5.5, and suppose
$d_{C^{1}}\left(\varphi_{j},\varphi\right)\rightarrow 0$ and
$d_{C^{1}}\left(g_{j},g\right)\rightarrow 0$ as $j\rightarrow\infty$. Then, it
holds that

${\displaystyle\prod_{t=1}^{T}p_{j}(y_{t}\mid
y^{t-1};\gamma)\overset{\mathit{l}_{\infty}}{\longrightarrow}\prod_{t=1}^{T}p(y_{t}\mid
y^{t-1};\gamma)}.$

The result shows that, as the researcher gets better approximations of the
policy function in a dynamic economic model, the computed likelihood converges
to the exact likelihood, even if the shock is unbounded. This result goes
beyond the result in Fernandez-Villaverde, Rubio-Ramirez, and Santos (2006)
and is particularly relevant to researchers using dynamic economic models with
unbounded shocks, such as normally distributed shocks– a standard
specification in the literature– as it guarantees, asymptotically, that the
likelihood function implied by the model is the correct object of interest.

## References

  * (1)
  * Dudley (2002) Dudley, R. M. (2002): “Real analysis and probability, volume 74 of Cambridge Studies in Advanced Mathematics,” .
  * Fernandez-Villaverde, Rubio-Ramirez, and Santos (2006) Fernandez-Villaverde, J., J. Rubio-Ramirez, and M. Santos (2006): “Convergence Properties of the Likelihood of Computed Dynamic Models,” _Econometrica_ , 74(1), 93–119.
  * Futia (1982) Futia, C. (1982): “Invariant distributions and the limiting behavior of Markovian economic models,” _Econometrica_ , 50(2), 377–408.
  * Kamihigashi (2007) Kamihigashi, T. (2007): “Stochastic optimal growth with bounded or unbounded utility and with bounded or unbounded shocks,” _Journal of Mathematical Economics_ , 43(3), 477–500.
  * Kamihigashi and Stachurski (2016) Kamihigashi, T., and J. Stachurski (2016): “Seeking ergodicity in dynamic economies,” _Journal of Economic Theory_ , 163, 900–924.
  * Munkres (2000) Munkres, J. R. (2000): _Topology_. Prentice Hall.
  * Nishimura and Stachurski (2005) Nishimura, K., and J. Stachurski (2005): “Stability of Optimal growth models: a new approach,” _Journal of Economic Theory_ , 122, 100–118.
  * Santos and Peralta-Alva (2005) Santos, M., and A. Peralta-Alva (2005): “Accuracy of simulations for stochastic dynamic models,” _Econometrica_ , 73(6), 1939–1976.
  * Santos (2004) Santos, M. S. (2004): “Simulation-based estimation of dynamic models with continuous equilibrium solutions,” _Journal of Mathematical Economics_ , 40(3-4), 465–491.
  * Smets and Wouters (2007) Smets, F., and R. Wouters (2007): “Shocks and frictions in US business cycles: A Bayesian DSGE approach,” _American economic review_ , 97(3), 586–606.
  * Stachurski (2002) Stachurski, J. (2002): “Stochastic Optimal growth with unbounded shock,” _Journal of Economic Theory_ , 106, 40–65.
  * Yosida and Kakutani (1941) Yosida, K., and S. Kakutani (1941): “Operator-theoretical treatment of Markoff’s process and mean ergodic theorem,” _Annals of Mathematics_ , pp. 188–228.

## Appendix

###### Proof of Lemma 3.1.

The “Only if” part is derived from $C^{\infty}\subset BL$. In this proof we
will show the “if” part. For any $f\in BL\left(S\right)$ and $\varepsilon>0$,
there exists $u\in C^{\infty}\left(S\right)$, such that $\sup_{s\in
S}\left|f\left(s\right)-u\left(s\right)\right|<\varepsilon$ by the mollifier
method. By triangular inequality, we have

$\displaystyle\sup_{f\in
BL\left(S\right)}\left|\int_{S}f\left(s\right)d\mu_{n}-\int_{S}f\left(s\right)d\mu\right|$
$\displaystyle\leq$ $\displaystyle\sup_{f\in
BL\left(S\right)}\left|\int_{S}f\left(s\right)d\mu_{n}-\int_{S}u\left(s\right)d\mu_{n}\right|+\sup_{u\in
C^{\infty}\left(S\right)}\left|\int_{S}u\left(s\right)d\mu_{n}-\int_{S}u\left(s\right)d\mu\right|$
$\displaystyle+\sup_{f\in
BL\left(S\right)}\left|\int_{S}u\left(s\right)d\mu-\int_{S}f\left(s\right)d\mu\right|.$

The first and third term of the right hand side of the inequality are smaller
than $\varepsilon$. The second term is given in the definition. ∎

###### Proof of Proposition 3.1.

The strong convergence of the sequence of operators is

$\displaystyle\sup_{s\in
S}\left\|Tf\left(s\right)-T_{j}f\left(s\right)\right\|$
$\displaystyle=\sup_{s\in S}\left|\int
f\left(\varphi\left(s,\varepsilon\right)\right)Q\left(s,d\varepsilon\right)-\int
f\left(\varphi_{j}\left(s,\varepsilon\right)\right)Q\left(s,d\varepsilon\right)\right|$
$\displaystyle=\sup_{s\in
S}\left|\int\left\\{f\left(\varphi\left(s,\varepsilon\right)\right)-f\left(\varphi_{j}\left(s,\varepsilon\right)\right)\right\\}Q\left(s,d\varepsilon\right)\right|$
$\displaystyle\rightarrow 0$

for all $f\in C^{2}\left(S\right)$. This simply means that

$\displaystyle\lim_{j\rightarrow\infty}\mathbb{E}\left[f\left(\varphi_{j}\left(s\right)\right)\right]$
$\displaystyle=\mathbb{E}\left[f\left(\varphi\left(s\right)\right)\right],$ $\
{}^{\forall}f\in C^{2}\left(S\right).$

Let $f$ belong to $\mathcal{A}$. Then, for any two $T$ and $T_{j}$, and
corresponding invariant measures $\mu$ and $\mu_{j}$, we have

$\displaystyle\left|\left\langle Tf,\hat{\mu}\right\rangle-\left\langle
T_{j}f,\mu_{j}\right\rangle\right|$ $\displaystyle\leq\left|\left\langle
Tf,\mu\right\rangle-\left\langle
Tf,\mu_{j}\right\rangle\right|+\left|\left\langle
Tf,\mu_{j}\right\rangle-\left\langle T_{j}f,\mu_{j}\right\rangle\right|.$ (8)

It follows from Proholov’s Theorem that $\left\\{\mu_{j}\right\\}$ has a
weakly convergent subsequence. Let $\left\\{\mu_{j_{k}}\right\\}$ be such a
subsequence, and let $\hat{\mu}$ be its limit. Then for the first term in eq.
(8), we have

$\displaystyle\left|\left\langle Tf,\hat{\mu}\right\rangle-\left\langle
Tf,\mu_{j_{k}}\right\rangle\right|$ $\displaystyle\rightarrow 0.$

Next, we show the equality: $\hat{\mu}=\mu_{0}$. We have

$\displaystyle\left|\left\langle f,\hat{\mu}\right\rangle-\left\langle
T_{0}f,\hat{\mu}\right\rangle\right|$ $\displaystyle\leq$
$\displaystyle\left|\left\langle f,\hat{\mu}\right\rangle-\left\langle
f,\mu_{j_{k}}\right\rangle\right|+\left|\left\langle
f,\mu_{j_{k}}\right\rangle-\left\langle
T_{0}f,\mu_{j_{k}}\right\rangle\right|$ $\displaystyle+\left|\left\langle
T_{0}f,\mu_{j_{k}}\right\rangle-\left\langle
T_{0}f,\hat{\mu}\right\rangle\right|.$

Since $Tf$ and $T_{0}f$ are continuous and $\left\\{\mu_{j_{k}}\right\\}$
convergence weakly to $\hat{\mu}$, the first and third terms on the right hand
side approaches zero as $j$ goes to infinity. For the second term, we have

$\displaystyle\left|\left\langle f,\mu_{j_{k}}\right\rangle-\left\langle
T_{0}f,\mu_{j_{k}}\right\rangle\right|$ $\displaystyle=\left|\left\langle
f,T_{j_{k}}^{*}\mu_{j_{k}}\right\rangle-\left\langle
T_{0}f,\mu_{j_{k}}\right\rangle\right|$ $\displaystyle=\left|\left\langle
T_{j_{k}}f,\mu_{j_{k}}\right\rangle-\left\langle
T_{0}f,\mu_{j_{k}}\right\rangle\right|$
$\displaystyle\leq\left\|T_{j_{k}}f-T_{0}f\right\|$

∎

###### Proof of Theorem 3.2.

The topology of weak convergence can be defined by the metric on the
probability measure space,

$\displaystyle d\left(\mu,\nu\right)$
$\displaystyle=\sup_{f\in\mathcal{A}}\left\\{\left|\int
f\left(s\right)\mu\left(ds\right)-\int
f\left(s\right)\nu\left(ds\right)\right|\right\\},$

where $\mathcal{A}$ is the space of Lipschitz functions on $S$, with constant
$L\leq 1$ and $-1\leq f\leq 1$. Then, we have

$\displaystyle\left\|Tf\left(s\right)-T_{j}f\left(s\right)\right\|$
$\displaystyle=\left|\int
f\left(\varphi\left(s,\varepsilon\right)\right)Q\left(d\varepsilon\right)-\int
f\left(\varphi_{j}\left(s,\varepsilon\right)\right)Q\left(d\varepsilon\right)\right|$
$\displaystyle=\left|\int
f\left(\varphi\left(s,\varepsilon\right)\right)-f\left(\varphi_{j}\left(s,\varepsilon\right)\right)Q\left(d\varepsilon\right)\right|.$

Since $f\in C^{\infty}$, there exists a constant $K$, such that

$\displaystyle\left|\int\left[f\left(\varphi\left(s,\varepsilon\right)\right)-f\left(\varphi_{j}\left(s,\varepsilon\right)\right)\right]Q\left(d\varepsilon\right)\right|$
$\displaystyle\leq Kd\left(\varphi,\varphi_{j}\right).$

∎

###### Proof of Proposition 4.1.

Denote $s_{n}\left(s_{0}\right)$ as

$\underbrace{\varphi\left(\varphi\left(\varphi\cdots\left(\varphi\left(s_{0},\varepsilon_{1}\right),\varepsilon_{2}\right)\right)\right)}_{\textrm{n
times}}.$

And denote $\hat{s}_{n}\left(s_{0}\right)$ as

$\underbrace{\hat{\varphi}\left(\hat{\varphi}\left(\hat{\varphi}\cdots\left(\hat{\varphi}\left(s_{0},\varepsilon_{1}\right),\varepsilon_{2}\right)\right)\right)}_{\textrm{n
times}}.$

Since $f$ is a Lipschitz function with constant $L$, we have

$\displaystyle\left|\mathbb{E}\left[f\left(s_{n}\left(s_{0}\right)\right)\right]-\mathbb{E}\left[f\left(\hat{s}_{n}\left(s_{0}\right)\right)\right]\right|$
$\displaystyle=$
$\displaystyle\left|\mathbb{E}\left[f\left(\varphi\left(s_{n-1}\left(s_{0}\right),\varepsilon_{n}\right)\right)\right]-\mathbb{E}\left[f\left(\hat{\varphi}\left(\hat{s}_{n-1}\left(s_{0}\right),\varepsilon_{n}\right)\right)\right]\right|$
$\displaystyle\leq$
$\displaystyle\left|\mathbb{E}\left[f\left(\varphi\left(s_{n-1}\left(s_{0}\right),\varepsilon_{n}\right)\right)\right]-\mathbb{E}\left[f\left(\varphi\left(\hat{s}_{n-1}\left(s_{0}\right),\varepsilon_{n}\right)\right)\right]\right|$
$\displaystyle+\left|\mathbb{E}\left[f\left(\varphi\left(s_{n-1}\left(s_{0}\right),\varepsilon_{n}\right)\right)\right]-\mathbb{E}\left[f\left(\hat{\varphi}\left(\hat{s}_{n-1}\left(s_{0}\right),\varepsilon_{n}\right)\right)\right]\right|$
$\displaystyle\leq$ $\displaystyle
L\left(1-\varepsilon\right)\mathbb{E}\left\|s_{n-1}\left(s_{0}\right)-\hat{s}_{n-1}\left(s_{0}\right)\right\|+Ld\left(\hat{\varphi},\varphi\right).$

By the same argument above, we have

$\displaystyle
L\left(1-\varepsilon\right)\mathbb{E}\left\|s_{n-1}\left(s_{0}\right)-\hat{s}_{n-1}\left(s_{0}\right)\right\|$
$\displaystyle\leq$ $\displaystyle
L\left(1-\varepsilon\right)^{2}\mathbb{E}\left\|s_{n-2}\left(s_{0}\right)-\hat{s}_{n-2}\left(s_{0}\right)\right\|+L\left(1-\varepsilon\right)^{2}d\left(\hat{\varphi},\varphi\right).$

Iterating this, we obtain

$\left|\mathbb{E}\left[f\left(s_{n}\left(s_{0}\right)\right)\right]-\mathbb{E}\left[f\left(\hat{s}_{n}\left(s_{0}\right)\right)\right]\right|\leq\frac{Ld\left(\hat{\varphi},\varphi\right)}{\varepsilon}.$

Integrating this by an invariant measure $\hat{\mu}^{*}$ of $\hat{\varphi}$,
we have

$\left|\int\mathbb{E}\left[f\left(s_{n}\left(s_{0}\right)\right)\right]\hat{\mu}^{*}\left(ds_{0}\right)-\int\mathbb{E}\left[f\left(\hat{s}_{n}\left(s_{0}\right)\right)\right]\hat{\mu}^{*}\left(ds_{0}\right)\right|\leq\frac{Ld\left(\hat{\varphi},\varphi\right)}{\varepsilon}.$

Note that the second term on the left-hand side is equal to $\int
f\left(s\right)\hat{\mu}^{*}\left(ds\right)$. From Theorem 4.1, for every
$s_{0}$, $\mathbb{E}\left[f\left(s_{n}\left(s_{0}\right)\right)\right]$
converges uniformly on $\beta\left(S\right)$ to $\int
f\left(s\right)\mu^{*}\left(ds\right)$. Finally, we obtain

$\left|\int f\left(s\right)\mu^{*}\left(ds\right)-\int
f\left(s\right)\hat{\mu}^{*}\left(ds\right)\right|\leq\frac{Ld\left(\hat{\varphi},\varphi\right)}{\varepsilon}.$

∎

######  Theorem .2.

(Change of variable formula). Let $g{:}\
\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$ be Lipschitz continuous, $n\leq m$.
Then, for each measurable function $p{:}\
\mathbb{R}^{n}\rightarrow\mathbb{R}$,

$\displaystyle\int_{\mathbb{R}^{n}}p\left(x\right)Jg\left(x\right)m\left(dx\right)$
$\displaystyle=\int_{\mathbb{R}^{n}}\left[\sum_{x\in
g^{-1}\left(y\right)}p\left(x\right)\right]m\left(dy\right).$

###### Proof of Proposition 5.1.

Let

$\displaystyle{\displaystyle\prod_{t=1}^{T}p_{j}(y_{t}\mid y^{t-1};\theta)}$
$\displaystyle=\int\left(\int\prod_{t=1}^{T}p(\varepsilon_{1}^{t};\theta)p_{j}(y_{t}\mid\varepsilon_{1}^{t},s_{0},y^{t-1};\theta)d\varepsilon_{1}^{t}\right)\mu_{j}^{*}(ds_{0};\theta),$
$\displaystyle{\displaystyle\prod_{t=1}^{T}\tilde{p}_{j}(y_{t}\mid
y^{t-1};\theta)}$
$\displaystyle=\int\left(\int\prod_{t=1}^{T}p(\varepsilon_{1}^{t};\theta)p(y_{t}\mid\varepsilon_{1}^{t},S_{0},y^{t-1};\theta)d\varepsilon_{1}^{t}\right)\mu_{j}^{*}(ds_{0};\theta),$
$\displaystyle{\displaystyle\prod_{t=1}^{T}p(y_{t}\mid y^{t-1};\theta)}$
$\displaystyle=\int\left(\int\prod_{t=1}^{T}p(\varepsilon_{1}^{t};\theta)p(y_{t}\mid\varepsilon_{1}^{t},s_{0},y^{t-1};\theta)d\varepsilon_{1}^{t}\right)\mu^{*}(ds_{0};\theta).$

Consider the convergence of

$\displaystyle
p_{j}\left(y_{t}\mid\varepsilon_{1}^{t},s_{0},y^{t-1};\theta\right)$
$\displaystyle=p(\eta_{j,t};\theta)p(\varepsilon_{j,2,t};\theta)\left|dy_{j}(\eta_{j,t},\varepsilon_{j,2,t};\theta)\right|,$

where

$\mid
dy_{j}(\mathit{v}_{t},\mathit{w}_{2,t};\gamma)\mid=det\left[\begin{array}[]{cc}\nabla\frac{\partial
g_{j}}{\partial v_{t}}&\nabla\frac{\partial g_{j}}{\partial
w_{2,t}}\end{array}\right].$

Since the change of variable formula eq. (.2) holds for
$Jg_{j}\left(\varepsilon,\eta;\theta\right)=\left|dy_{j}(\mathit{v}_{j,t},\mathit{w}_{j,2,t};\theta)\right|$,
we have

$\int p_{j}\left(y_{t}\mid
W_{1}^{t},S_{0},y^{t-1};\theta\right)Q\left(d\varepsilon,d\eta\right)=\int
p\left(\varepsilon;\theta\right)p\left(\eta;\theta\right)Jg_{j}\left(\varepsilon,\eta;\theta\right)Q\left(d\varepsilon,d\eta\right).$

By $d_{C^{1}}\left(\varphi_{j},\varphi\right)\rightarrow 0$ and
$d_{C^{1}}\left(g_{j},g\right)\rightarrow 0$, we have

$\int
p\left(\varepsilon;\theta\right)p\left(\eta;\theta\right)Jg_{j}\left(\varepsilon,\eta;\theta\right)Q\left(d\varepsilon,d\eta\right)\rightarrow\int
p\left(\varepsilon;\theta\right)p\left(\eta;\theta\right)Jg\left(\varepsilon,\eta;\theta\right)Q\left(d\varepsilon,d\eta\right).$

Then, we have

$\prod_{t=1}^{T}p_{j}(y_{t}\mid
y^{t-1};\theta)\rightarrow{\displaystyle\prod_{t=1}^{T}\tilde{p}_{j}(y_{t}\mid
y^{t-1};\theta)}.$

By Lemma 1 of Fernandez-Villaverde, Rubio-Ramirez, and Santos (2006),
$\left(\int\prod_{t=1}^{T}p(W_{1}^{t};\theta)p(y_{t}\mid
W_{1}^{t},S_{0},y^{t-1};\theta)dW_{1}^{t}\right)$ is a continuous function.
Further, given Proposition 1 of Fernandez-Villaverde, Rubio-Ramirez, and
Santos (2006), we have
$\mu_{j}^{*}(dS_{0};\theta)\rightarrow\mu^{*}(dS_{0};\theta)$. Therefore,

$\displaystyle{\displaystyle\prod_{t=1}^{T}\tilde{p}_{j}(y_{t}\mid
y^{t-1};\theta)}$
$\displaystyle\rightarrow{\displaystyle\prod_{t=1}^{T}p(y_{t}\mid
y^{t-1};\theta)}.$

∎


# Optimizing carbon tax for decentralized electricity markets using an agent-
based model

Alexander J. M. Kell, A. Stephen McGough, Matthew Forshaw  School of
ComputingNewcastle UniversityNewcastle upon TyneUK a.kell2, stephen.mcgough,
<EMAIL_ADDRESS>
(2020)

###### Abstract.

Averting the effects of anthropogenic climate change requires a transition
from fossil fuels to low-carbon technology. A way to achieve this is to
decarbonize the electricity grid. However, further efforts must be made in
other fields such as transport and heating for full decarbonization. This
would reduce carbon emissions due to electricity generation, and also help to
decarbonize other sources such as automotive and heating by enabling a low-
carbon alternative. Carbon taxes have been shown to be an efficient way to aid
in this transition.

In this paper, we demonstrate how to to find optimal carbon tax policies
through a genetic algorithm approach, using the electricity market agent-based
model ElecSim. To achieve this, we use the NSGA-II genetic algorithm to
minimize average electricity price and relative carbon intensity of the
electricity mix. We demonstrate that it is possible to find a range of carbon
taxes to suit differing objectives.

Our results show that we are able to minimize electricity cost to below
£10/MWh as well as carbon intensity to zero in every case. In terms of the
optimal carbon tax strategy, we found that an increasing strategy between 2020
and 2035 was preferable. Each of the Pareto-front optimal tax strategies are
at least above £81/tCO2 for every year. The mean carbon tax strategy was
£240/tCO2.

Energy markets, policy, carbon tax, genetic algorithm, optimization, digital
twin, agent-based models, electricity market model, climate change

††journalyear: 2020††copyright: acmlicensed††conference: The Eleventh ACM
International Conference on Future Energy Systems; June 22–26, 2020; Virtual
Event, Australia††booktitle: The Eleventh ACM International Conference on
Future Energy Systems (e-Energy’20), June 22–26, 2020, Virtual Event,
Australia††price: 15.00††doi: 10.1145/3396851.3402369††isbn:
978-1-4503-8009-6/20/06††ccs: Computing methodologies Agent / discrete
models††ccs: Computing methodologies Genetic algorithms

##  1\. Introduction

Computer simulation allows practitioners to model real-world systems using
software. These simulations allow for ‘what-if ’ analyses which can provide an
indication as to how a system may behave under specific policies, environments
and assumptions. These simulations become particularly important in real
systems which have high costs, impacts or risks associated with them. A
‘digital twin’ is a concept which has emerged in recent years. This is defined
as a simulation of a specific instance of a system. This digital twin can then
be used to learn of optimization techniques that can be applied to the real
system. This foregoes the necessity of experimenting with the real system,
avoiding potential adverse side effects.

An electricity market is an example of a complex system which can be modelled
using a digital twin. Disruptions to electricity supply, a substantial
increase in the cost of electricity or unrestrained carbon emissions have the
potential to destabilize economies (Kaseke and Hosking, 2013; Masson-Delmotte
et al., 2018). It is for reasons such as these that electricity market models
are used to test hypotheses, develop strategies and gain an understanding of
underlying dynamics to prevent undesirable consequences (Jebaraj and Iniyan,
2006).

In this paper, we use the electricity market agent-based model ElecSim to find
an optimum carbon tax policy (Kell et al., 2019). Specifically, we use a
genetic algorithm to find a carbon tax policy to reduce both average
electricity price and the relative carbon density by 2035 for the UK
electricity market.

Carbon taxes have been shown to quickly lower emissions and lower the costs to
the public (Wittneben, 2009). Carbon taxes are able to send clear price
signals, as opposed to a cap-and-trade scheme, such as the EU Emissions
Trading System, which has shown to be unstable (Wittneben, 2009).

In this paper, we use the reference scenario projected by the UK Government’s
Department for Business & Industrial Strategy (BEIS) with model parameters
calibrated by Kell et al. (Department for Business Energy & Industrial
Strategy, 2019; Kell et al., 2020). This reference scenario projects energy
and emissions until 2035. We undertake various carbon tax policy interventions
to see how we could reduce relative carbon density whilst at the same time,
reduce the average electricity price.

The parameter space we optimize over is the carbon tax price over a 17 year
period from 2018 to 2035. The carbon price is used to influence the objectives
of average electricity price and relative carbon intensity in 2035. Grid and
random search are approaches which trial parameters at evenly distributed
spaces and random spaces respectively. These approaches are often inefficient,
however, and require an increased number of simulations due to their static
nature. Genetic Algorithms, in contrast, use an evolutionary computing
approach to find global optimal solutions faster.

In order to optimize over two potentially competing objectives, i.e. average
electricity price and relative carbon intensity, we use the Non-Dominated
Sorting Genetic Algorithm II (NSGA-II) (Deb et al., 2000). The NSGA-II
algorithm can approximate a Pareto frontier (Pareto and Schwier, 1927;
Stadler, 1979). A Pareto frontier is a curve in which there is no solution
which is better than another along the curve for different sets of parameters.
In this context, better means that a solution is closer to the optimal for a
particular combination of objectives.

We find that the rewards of average electricity price and relative carbon
intensity are not mutually exclusive. That is, it is possible to have both a
lower average electricity price and a lower relative carbon price. This is due
to the low short-run marginal cost of renewable technology, which has been
shown to lower electricity prices (O’Mahoney and Denny, 2011).

The main contribution of this paper is to explore carbon tax strategies using
genetic algorithms for multi-objective optimization.

The following sections are set out as follows. Section 2 covers examples of
optimisations using genetic algorithms and different carbon strategies.
Section 3 details the optimization techniques applied. Section 4 explores the
electricity market model used. We present our results in Section 5, and
conclude in Section 6.

##  2\. Literature Review

Multi-objective optimization problems are commonplace. In this section, we
review multiple applications that have used multi-objective optimization, as
well as explore the literature which focus on finding optimal carbon tax
strategies.

###  2.1. Examples of Optimization

Similar to our work, Ascione et al. use the NSGA-II algorithm to generate a
Pareto front to optimize for two objectives: operating cost for space
conditioning and thermal comfort (Ascione et al., 2016). The aim of their
paper is to optimize the hourly set point temperatures with a day-ahead
planning horizon. A Pareto front is generated, which allows a user to select a
solution according to their comfort needs and economic constraints. This work
showed a reduction in operating costs by up to 56% as well as improved thermal
comfort.

Gorzałczany et al. also apply the NSGA-II algorithm. However, they apply it to
the credit classification problem (Gorzalczany and Rudziński, 2016). The
objectives optimized over were accuracy and interpretability when making
financial decisions such as credit scoring and bankruptcy prediction. This
technique was able to significantly outperform the alternative methods in
terms of interpretability while remaining competitive or superior in terms of
the accuracy and speed of decision making in comparison with the existing
classification methods.

###  2.2. Carbon Tax Strategies

In this section, we explore different strategies employed in the literature to
analyze the benefits and consequences of a carbon tax. To the best of our
knowledge, we are the first to employ a multi-objective optimization algorithm
to minimize average electricity price and relative carbon density.

Levin et al. use an optimization model to analyze market and investment
impacts of several incentive mechanisms to support investment in renewable
energy and carbon emission reductions (Levin et al., 2019). Carbon tax was
found to be the most cost-efficient method of reducing emissions.

Zhou et al. construct a social welfare model based on a Stackelberg game (Zhou
et al., 2019). The differences and similarities between a flat carbon tax and
an increasing block tariff carbon tax are analyzed using a numerical
simulation. This work shows that an increasing block tariff carbon tax policy
can significantly reduce tax burdens for manufacturers and encourage low-
carbon production. In contrast to Zhou et al. we trial multiple different
carbon tax strategies using a machine learning approach.

Li et al. use a hierarchical carbon market scheduling model to reduce carbon
emissions (Li et al., 2017). Multi-objective optimization was applied to
discover optimal behaviours for policymakers, customers and generators to
minimize the costs incurred by these actors. Our work, however, focuses on the
different strategies of carbon tax as opposed to optimal actor behaviour.

##  3\. Optimization methods

Multi-objective optimization allows practitioners to overcome the problems
with optimizing multiple objectives with classical optimization techniques.
Multi-objective optimization algorithms are able to generate Pareto-optimal
solutions as opposed to converting the multiple objectives into a single-
objective problem. A single-objective problem assumes that there is only a
single optimum, and that other combinations are inferior. This may not be the
case, as different solutions are superior for a different set of
circumstances. A Pareto frontier is made up of many Pareto-optimal solutions
which can be displayed graphically. A user is then able to choose between
various solutions and trade-offs according to their wishes.

The NSGA-II algorithm, a multi-objective genetic optimization algorithm, is
able to generate a Pareto frontier in a single optimization run.

In the following sub-sections, we detail the standard genetic algorithm
followed by the NSGA-II algorithm.

###  3.1. Genetic Algorithms

Genetic Algorithms (GAs) (Holland, 1975) are a class of evolutionary algorithm
which can be used for optimization.

As shown by Algorithm 1, initially, a population of parameters $P_{0}$ are
generated to be trialled in the simulation. $P_{0}$ is then evaluated for
fitness, by running an instance of the simulation with the respective
parameters as inputs. Where fitness in this case is the reward that is to be
optimised. Next, a subset of individuals from $P_{0}$ are chosen for mating,
$C_{t+1}\subset P_{t}$. This subset of individuals are selected proportional
to their fitness. For mating with the subset $C_{t+1}$, the ‘fitter’
individuals have a higher chance of reproducing to create the offspring group
$C^{\prime}_{t+1}$. The individuals of $C^{\prime}_{t+1}$ have characteristics
dependent on the genetic operators, crossover and mutation (Mitchell, 1998).
These genetic operators are an implementation decision (Back et al., 2009).
The new population $P_{t+1}$ is then created by merging individuals from
$C^{\prime}_{t+1}$ and $P_{t}$. See Algorithm 1 for detailed pseudocode.

1:$t=0$

2:initialize $P_{t}$

3:evaluate structures in $P_{t}$

4:while termination condition not satisfied do

5: $t=t+1$

6: select reproduction $C_{t}$ from $P_{t-1}$

7: recombine and mutate structures in $C_{t}$ forming $C^{\prime}_{t}$

8: evaluate structures in $C^{\prime}_{t}$

9: select each individual for $P_{t}$ from $C^{\prime}_{t}$ or $P_{t-1}$

10:end while

Algorithm 1 Genetic algorithm (Back et al., 2009)

###  3.2. NSGA-II

NSGA-II is efficient for multi-objective optimization on a number of benchmark
problems and finds a better spread of solutions than Pareto Archived Evolution
Strategy (PAES) (Knowles and Corne, 1999) and Strength Pareto EA (SPEA)
(Zitzler and Thiele, 2006) when approximating the true Pareto-optimal front
(Deb et al., 2000).

The majority of multi-objective optimization algorithms use the concept of
_domination_ during population selection (Burke and Graham, 2014). A non-
dominated algorithm, however, seeks to achieve the Pareto-optimal solution.
This is where no single solution should dominate another. An individual
solution $\mathbf{x}^{1}$ is said to dominate another $\mathbf{x}^{2}$, if and
only if there is no objective of $\mathbf{x}^{1}$ that is worse than objective
of $\mathbf{x}^{2}$ and at least one objective of $\mathbf{x}^{1}$ is better
than the same objective of $\mathbf{x}^{2}$ (Bao et al., 2017).

Non-domination sorting is the process of finding a set of solutions which do
not dominate each other and make up the Pareto front. See Figure 1a for a
visual representation, where $f_{1}$ and $f_{2}$ are two objectives to
minimize and L1, L2 and L3 are dominated layers.

In this section, we define the processes used by the NSGA-II algorithm to
determine which solutions to keep:

####  3.2.1. Non-dominated sorting

We assume that there are $M$ objective functions to minimise, and that ${\bf
x^{1}}=\\{x_{j}^{1}\\}$ and $\bf x^{2}$ are two solutions.
$x_{j}^{1}<x_{j}^{2}$ implies solution $\bf x^{1}$ is better than solution
$\bf x^{2}$ on objective $j$. A solution $\mathbf{x}^{1}$ is said to dominate
the solution $\mathbf{x}^{2}$ if the following conditions are true:

  1. (1)

The solution $\mathbf{x}^{1}$ is no worse than $\mathbf{x}^{2}$ in every
objective. I.e. $x^{1}_{j}\leq x^{2}_{j}\;\;\forall j\in\\{1,2,\ldots,M\\}$.

  2. (2)

The solution $\mathbf{x}^{1}$ is better than $\mathbf{x}^{2}$ in at least one
objective. I.e. $\exists\
{j}\in\\{1,2,\ldots,M\\}\;\;s.t.\;\;x^{1}_{j}<x^{2}_{j}$.

Figure 1. a) Schematic of non-dominated sorting with solution layering b)
Schematic of the NSGA-II procedure

Next, each of the solutions are ranked according to their level of non-
domination. An example of this ranking is shown in Figure 1a. Here, $f_{1}$
and $f_{2}$ are the objectives to be minimised. The Pareto front is the first
front. All of the solutions in the Pareto front are not dominated by any other
solution. The solutions in layer 1, L1, are dominated only by those in the
Pareto front, and are non-dominated by those in L2 and L3.

The solutions are then ranked according to their layer. For example, the
solutions in the Pareto front are given a fitness rank ($i_{rank}$) of 1,
solutions in L1 have an $i_{rank}$ of 2.

####  3.2.2. Density Estimation

($i_{distance}$) is calculated for each solution. This is the average distance
between the two closest points to the solution in question.

####  3.2.3. Crowded comparison operator

($\prec_{n}$) is used to ensure that the final frontier is an evenly spread
out Pareto-optimal front. This is achieved by using the two attributes:
$(i_{rank})$ and$(i_{distance})$. The partial order is then defined as:
$i\prec_{n}j$ if $(i_{rank}<j_{rank})$ or $((i_{rank}=j_{rank})$ and
$(i_{distance}>j_{distance}))$ (Deb et al., 2000).

This order prefers solutions with a lower rank $i_{rank}$. For solutions with
the same rank, the solution in the less dense area is preferred.

####  3.2.4. Main loop

Similarly to the standard GA, a population $P_{0}$ is created with random
parameters. The solutions of $P_{0}$ are then sorted according to non-
domination. The child population $C^{\prime}_{1}$ of size $N$ is then created
by binary tournament selection, recombination and mutation operators. In this
case, tournament selection is the process of evaluating and comparing the
fitness of various individuals within a population. In binary tournament
selection, two individuals are chosen at random, the fitnesses are evaluated,
and the individual with the better solution is selected (Abd Rahman et al.,
2016).

1:$R_{t}=P_{t}\cup C^{\prime}_{t}$ combine parent and child population

2:$\mathcal{F}=$ fast-non-dominated-sort $(R_{t})$ where
$\mathcal{F}=(\mathcal{F}_{1},\mathcal{F}_{2},\ldots)$

3:$P_{t+1}=\emptyset$

4:while $\ \textbf{do}\left|P_{t+1}<N\right|$

5: Calculate the crowding distance of $(\mathcal{F}_{i})$)

6: $P_{t+1}=P_{t+1}\cup\mathcal{F}_{i}$

7:end while

8:Sort($P_{t+1},\prec_{n}$) sort in descending order using $\prec_{n}$

9:$P_{t+1}=P_{t+1}[0:N]$ select the first $N$ elements of $P_{t+1}$

10:$Q_{t+1}=$ make-new-population$(P_{t+1})$ using selection, crossover and
mutation to create the new population $Q_{t+1}$

11:$t=t+1$

Algorithm 2 NSGA-II main loop (Deb et al., 2000)

Next, a new combined population is formed $R_{t}=P_{t}\cup C^{\prime}_{t}$.
$R_{t}$ has a size of $2N$. $R_{t}$ is then sorted according to non-
domination. A new population is then formed $P_{t+1}$. Solutions are added
from the sorted $R_{t}$ in order of non-domination. Solutions are added until
the size of $P_{t+1}$ exceeds $N$. The solutions from the last layer are
prioritised based on having the largest crowding distance (Deb et al., 2000).

This process is shown in Figure 1b, which is repeated until the termination
condition is met. A termination condition could be: no significant improvement
over $X$ iterations or a specified number of iterations have been performed.
The full procedure is detailed formally by Algorithm 2.

###  3.3. Carbon Optimization Application

In this section, we describe how the genetic algorithm is applied in our
carbon optimization case. We use multi-objective optimization to find a
solution which has both a low carbon emission and low average electricity
price. The parameters that we adjust is the carbon tax between the years 2018
and 2035.

The mating steps work by, initially, taking the sets of carbon prices over the
17 year period (2018 to 2035) which have the best rewards (lowest relative
carbon emissions and average electricity price). These carbon prices are then
mated with a probability of 90%, creating child carbon prices. The children
are mutated with a probability of 5%. Therefore, 5% of children have a carbon
price which is not inherited from the parents. Over time, the mutations and
inherited properties tend to a population with more desirable rewards.

##  4\. Experimental Setup

###  4.1. Simulation Environment

In order to evaluate the different carbon strategies, we used the model
developed by Kell et al., ElecSim (Kell et al., 2019, 2020). ElecSim is an
agent-based model which mimics the behaviour of decentralized electricity
markets. For this paper, we parametrized the model to data for the UK in 2018
to act as a digital twin of the UK electricity market. This includes the power
plants in operation in 2018, and the funds available to their respective
companies (for Business Energy and Strategy, 2019; Department for Business,
2019). ElecSim is validated by being instantiated by data from 2013 and
projected forward to 2018, with a mean absolute scaled error (MASE) below or
equal to 0.701 for all generator types (Kell et al., 2020).

Six fundamental sections make up ElecSim: 1) power plant data; 2) scenario
data; 3) the time-steps of the algorithm; 4) the power exchange; 5) the
investment algorithm and 6) the generation companies (GenCos) as agents.
ElecSim uses a subset of representative days of electricity demand, solar
irradiance and wind speed to approximate a full year. In this context,
representative days are a subset of days which, when scaled up, can adequately
represent a year. Figure 2 details how these components interact.

Specifically, the configuration file details the scenario which can be set by
the user. This includes electricity demand, carbon price and fuel prices. The
data sources parametrize the digital twin to a particular country, including
information such as wind capacity and power plants in operation. Generation
Companies own and invest in power plants. These power plants are then matched
to electricity demand using a spot market.

Figure 2. System overview of ElecSim (Kell et al., 2019).

The market runs in merit-order dispatch and bids are made by the power plant’s
short-run marginal cost (SRMC). Investment in power plants are based upon a
net present value (NPV) calculation. Decisions are made every year, where new
generators enter the market. NPV is able to evaluate and compare investments
with cash flows spread over many years. This is shown formally in Equation 1,
where $t$ is the year of the cash flow, $i$ is the discount rate, $N$ is the
total number of years, or lifetime of the power plant, and $R_{t}$ is the net
cash flow of the year $t$:

(1) $NPV(t,N)=\sum_{t=0}^{N}\frac{R_{t}}{(1+t)^{t}}.$

$R_{t}$ takes into account the expected carbon price ten years into the
future. Generator companies use a linear regression based on historical data
to project the carbon price. With an increasing carbon price, the $R_{t}$
reduces for high carbon-emitting power plants, therefore making low carbon
options more desirable. As demand in ElecSim is exogenous, demand does not
change with respect to price and is assumed to be perfectly price inelastic.

The yearly income for each power plant is estimated for each generation
company by running a merit-order dispatch electricity market ten years into
the future. However, the expected cost of electricity ten years into the
future is uncertain. We, therefore, use the reference scenario projected by
BEIS and use the predicted costs of electricity calibrated by Kell et al
(Department for Business Energy & Industrial Strategy, 2019; Kell et al.,
2020). The agents predict the future carbon price by using a linear regression
model.

###  4.2. Optimization

In this section, we detail the optimization approach taken. We modify the
carbon tax each year, as we believe this is the most likely process taken by
governments, giving generator companies and consumers the ability to
understand market conditions during each year.

####  4.2.1. Non-parametric carbon policy

The optimization approach has two stages. First, we initialize the population
of the NSGA-II algorithm $P_{0}$ with 18 attributes. These correspond to a
separate carbon tax for each year, shown by Equation 2:

(2) $P_{0}=\\{a_{1},a_{2},\ldots,a_{18}\\},0\leq a_{y}\leq 250,$

where $P_{0}$ is the first population, $a_{y}$ is the attribute or carbon
price in year $y$ and $a_{1}$ is the carbon price in year 1, $a_{2}$ the
carbon price in year two and so forth. The constraints of the algorithm are
that each of the carbon prices are bound between the values of £$0$ and
£$250$. This provides the optimization algorithm with the highest degree of
freedom. The value £$250$ was chosen due to the relative costs of electricity,
where £$250$ would be the upper bound for the cost of electricity. This high
degree of freedom enables a high number of strategies to be trialled due to
its non-parametric nature. This, however, comes with a large search space
requiring a large number of iterations.

####  4.2.2. Linear carbon policy

To reduce the search space for the carbon strategy, we also trial a linear
carbon strategy, of the form:

(3) $p_{c}=a_{1}y_{t}+a_{2},-14\leq a_{1}\leq 14,0\leq a_{2}\leq 250,$

where $p_{c}$ is the carbon price, $y_{t}$ is the year, $a_{1}$ is the
gradient or first attribute and $a_{2}$ is the intercept or second attribute.
The constraints of the optimisation problem are that $a_{1}$ is bound by $-14$
and $14$, and $a_{2}$ by 0 and 250. These bounds are chosen to ensure that the
carbon price does not exceed ${\sim}$£500 in the year 18 (2035) and is greater
than about -£250, as well as ensuring that the carbon tax in the first year is
greater than £0 but smaller than £250. The bounds for $a_{1}$ was chosen to
make the mathematics simpler, whilst remaining in range.

##  5\. results

In this section, we explore the results of the optimizations, the optimum
carbon strategies and the resultant electricity mixes.

###  5.1. Non-parametric carbon policy

Figure 4 displays the development of the genetic algorithm against the
rewards, relative carbon density and average electricity price. Darker colours
display higher generation numbers. The first generation shows a widespread in
relative carbon density and average electricity price. However, over
successive generations, the solutions converge to a relative carbon density of
0 and an average electricity price under £10MW/h.

Strikingly, the rewards of relative carbon density and average electricity
price are not mutually destructive. This could be due to the low short-run
marginal cost of renewable energy which reduces both electricity prices and
carbon emissions (O’Mahoney and Denny, 2011).

To understand the optimum carbon strategies, we visualized the parameters that
produced the lowest average electricity prices in Figure 4. Specifically, we
filtered for electricity prices under £5/MWh and displayed the results using a
heat map. The darker colours represent a higher density of points.

Figure 4 displays a general trend, where carbon tax starts at ${\sim}$£100
until the year 2030, where it increases to ${\sim}$£200 by 2035. This may be
due to the fact that a lower initial carbon tax of ${\sim}$£100 encourages
investment in low-carbon technologies before the higher rate of ${\sim}$£200
comes into force. This higher rate of carbon tax would allow GenCos to
outcompete higher carbon-emitting generators over the lifetime of the plants.

Figure 3. Development of genetic algorithm rewards for non-parametric carbon
tax policy results in 2035.

Figure 4. Density plot of points with an average electricity price ¡£5/MWh for
non-parametric carbon tax policy results in 2035.

###  5.2. Linear carbon policy

Figure 5 displays the development of the genetic algorithm against the
rewards: relative carbon density and average electricity price. Similarly to
the non-parametric carbon policy shown in Figure 4, the first generation shows
a wide spread of results. However, the spread is smaller than that of the
linear carbon policy. This may be due to the fact that it is easier for the
GenCos to predict the carbon policy, which increases confidence in the NPV
calculations. The linear carbon policy also converges to a relative carbon
density of 0, and an average electricity price smaller than £10MW/h.

Figure 6 compares the distributions of average electricity price for both
techniques. Both methods show improvements as the number of generations of the
genetic algorithm increase. The linear policy, however, is able to more
quickly converge to a low average electricity price, with a mode of
${\sim}$£5.4MW/h. The non-parametric policy has a number of poorer performing
parameters, and Generation Number 4 has a bimodal distribution, with a mode of
${\sim}$£6.3MW/h.

Figure 5. Development of genetic algorithm rewards in 2035 for linear carbon
strategy. Figure 6. Density plot of average electricity price in 2035 over
generation number of genetic algorithm for both linear and non-parametric
policy.

Figure 7(a) displays the linear carbon policies which had an average
electricity price under £4.5MW/h. There is no single ‘optimum’ carbon policy;
a range of policies are able to achieve low carbon and a low average
electricity price.

We explore the electricity mix generated of three different strategies shown
in Figure 7(b). We selected the highest, lowest, and the lowest flat carbon
strategy to show a range of possible strategies.

(a) All carbon policies.

(b) Highlighted policies.

Figure 7. Linear carbon policies under £4.5MW/h visualised.

Figure 8 displays the generated electricity mixes for each of the selected
strategies. To generate these images, we ran 80 scenarios to capture the
variability between scenarios.

Whilst there does not seem to be a significant difference between scenarios,
with solar providing ${\sim}60\%$ of the electricity mix by 2035, there is an
observable difference with the other generator types.

The ‘highest’ carbon strategy exhibits a higher uptake in nuclear, possibly
due to the fact that nuclear becomes more competitive when compared to coal or
gas. The ‘lowest’ carbon strategy shows a higher uptake in Combined Cycle Gas
Turbines (CCGT) during the years of 2026 to 2031 as it outcompetes nuclear.
The ‘flat’ carbon policy shows a higher percentage of solar energy than any of
the other scenarios, albeit with a lower percentage of nuclear. Onshore wind
is shown to be consistent for these scenarios.

Figure 8. Electricity mixes under selected linear carbon policies.

##  6\. Conclusion

In this paper, we have demonstrated that it is possible to use the genetic
algorithm technique NSGA-II to optimize carbon tax policy using an electricity
market agent-based model.

We trialled a non-parametric carbon policy by allowing the genetic algorithm
to optimize a carbon price for each year. These results showed us that a
linear carbon tax might be appropriate. We then used a linear model as a
carbon tax policy to reduce the total number of parameters for the genetic
algorithm to optimize.

We were able to show that a range of linear carbon taxes were able to achieve
both low average electricity price and a relative carbon intensity of zero in
2035. By exploring three different carbon tax policies, we saw that
${\sim}$60% of electricity consumption in the UK would be provided by solar.
The difference between these ‘optimal’ carbon tax policies was largely shown
by competition between CCGT, coal and nuclear.

This was largely due to the low short-run marginal cost of solar and nuclear
energy, which means that they are often dispatched ahead of the fossil-fuel
based generators. CCGT and coal, however, are useful for filling demand when
there is low solar irradiance.

In future work, we would like to try additional scenarios with varying future
generation costs and calculate a sensitivity analysis to carbon taxes. In
addition to this, we would like to model the uncertain reactions by consumers
and generation companies with regards to carbon taxes. The linear carbon tax
approach is an introductory approach which can be expanded upon.

###### Acknowledgements.

This work was supported by the Engineering and Physical Sciences Research
Council, Centre for Doctoral Training in Cloud Computing for Big Data [grant
number EP/L015358/1].

## References

  * (1)
  * Abd Rahman et al. (2016) Rosshairy Abd Rahman, Razamin Ramli, Zainoddin Jamari, and Ku Ruhana Ku-Mahamud. 2016.  Evolutionary Algorithm with Roulette-Tournament Selection for Solving Aquaculture Diet Formulation.  _Mathematical Problems in Engineering_ 2016 (2016), 1–10.  https://doi.org/10.1155/2016/3672758
  * Ascione et al. (2016) Fabrizio Ascione, Nicola Bianco, Claudio De Stasio, Gerardo Maria Mauro, and Giuseppe Peter Vanoli. 2016.  Simulation-based model predictive control by the multi-objective optimization of building energy performance and thermal comfort.  _Energy and Buildings_ 111 (2016), 131–144.  https://doi.org/10.1016/j.enbuild.2015.11.033
  * Back et al. (2009) Thomas Back, David B Fogel, and Zbigniew Michalewicz. 2009\.  Evolutionary Computation 1 Basic Algorithms and Operators.  _Comprehensive Chemometrics_ (2009).
  * Bao et al. (2017) Chunteng Bao, Lihong Xu, Erik D. Goodman, and Leilei Cao. 2017\.  A novel non-dominated sorting algorithm for evolutionary multi-objective optimization.  _Journal of Computational Science_ 23 (2017), 31–43.  https://doi.org/10.1016/j.jocs.2017.09.015
  * Burke and Graham (2014) Edmund K. Burke and Kendall Graham. 2014.  _Search methodologies: Introductory tutorials in optimization and decision support techniques, second edition_.  1–716 pages.  https://doi.org/10.1007/978-1-4614-6940-7
  * Deb et al. (2000) Kalyanmoy Deb, Samir Agrawal, Amrit Pratap, and T Meyarivan. 2000.  A fast elitist non-dominated sorting genetic algorithm for multi-objective optimisation: NSGA-II.  _CEUR Workshop Proceedings_ 1133 (2000), 850–857.
  * Department for Business (2019) Energy & Industrial Strategy Department for Business. 2019.  Companies House - GOV.UK.  _UK Government_ (2019).
  * Department for Business Energy & Industrial Strategy (2019) Department for Business Energy & Industrial Strategy. 2019.  Updated energy and emissions projections 2018.  _The Energy White Paper_ April (2019).
  * for Business Energy and Strategy (2019) Department for Business Energy and Industrial Strategy. 2019.  Power stations in the United Kingdom, May 2019.  _Digest of United Kingdom Energy Statistics (DUKES)_ (2019).
  * Gorzalczany and Rudziński (2016) Marian B. Gorzalczany and Filip Rudziński. 2016.  A multi-objective genetic optimization for fast, fuzzy rule-based credit classification with balanced accuracy and interpretability.  _Applied Soft Computing Journal_ 40 (2016), 206–220.  https://doi.org/10.1016/j.asoc.2015.11.037
  * Holland (1975) John H. Holland. 1975\.  _Search methodologies: Introductory tutorials in optimization and decision support techniques, 2nd ed_.
  * Jebaraj and Iniyan (2006) S. Jebaraj and S. Iniyan. 2006.  A review of energy models.  _Renewable and Sustainable Energy Reviews_ 10, 4 (2006), 281–311.  https://doi.org/10.1016/j.rser.2004.09.004
  * Kaseke and Hosking (2013) Nyasha Kaseke and Stephen G Hosking. 2013.  Sub-Saharan Africa Electricity Supply Inadequacy: Implications.  29, 2 (2013), 113–132.  https://doi.org/10.1353/eas.2013.0009
  * Kell et al. (2019) Alexander Kell, Matthew Forshaw, and A Stephen Mcgough. 2019\.  ElecSim : Monte-Carlo Open-Source Agent-Based Model to Inform Policy for Long-Term Electricity Planning.  _The Tenth ACM International Conference on Future Energy Systems (ACM e-Energy)_ (2019), 556–565.
  * Kell et al. (2020) Alexander J. M. Kell, Matthew Forshaw, and A. Stephen McGough. 2020.  Long-Term Electricity Market Agent Based Model Validation using Genetic Algorithm based Optimization.  _The Eleventh ACM International Conference on Future Energy Systems (e-Energy’20)_ (2020).
  * Knowles and Corne (1999) Joshua Knowles and David Corne. 1999.  The Pareto archived evolution strategy: A new baseline algorithm for Pareto multiobjective optimisation.  _Proceedings of the 1999 Congress on Evolutionary Computation, CEC 1999_ 1 (1999), 98–105.  https://doi.org/10.1109/CEC.1999.781913
  * Levin et al. (2019) Todd Levin, Jonghwan Kwon, and Audun Botterud. 2019\.  The long-term impacts of carbon and variable renewable energy policies on electricity markets.  _Energy Policy_ 131, February (2019), 53–71.  https://doi.org/10.1016/j.enpol.2019.02.070
  * Li et al. (2017) Dan Li, Weiqi Hua, Hongjian Sun, and Wei Yu Chiu. 2017\.  Multiobjective optimization for carbon market scheduling based on behavior learning.  _Energy Procedia_ 142 (2017), 2089–2094.  https://doi.org/10.1016/j.egypro.2017.12.581
  * Masson-Delmotte et al. (2018) V Masson-Delmotte, P Zhai, H.O Pörtner, D Roberts, J Skea, P R Shukla, A Pirani, W Moufouma-Okia, C Péan, R Pidcock, S Connors, J B Matthews, Y Chen, X Zhou, M I Gomis, E Lonnoy, T Maycock, M Tignor, and T Waterfield. 2018.  _IPCC Special Report 1.5 - Summary for Policymakers_.  IPCC.
  * Mitchell (1998) Melanie Mitchell. 1998\.  _An introduction to genetic algorithms_.  MIT press.
  * O’Mahoney and Denny (2011) Amy O’Mahoney and Eleanor Denny. 2011.  The Merit Order Effect of Wind Generation on the Irish Electricity Market.  _Usaee 2011_ 56043 (2011), 1–12.  https://doi.org/10.1227/01.NEU.0000349921.14519.2A arXiv:suresh govindarajan
  * Pareto and Schwier (1927) Vilfredo Pareto and Ann S (Translator) Schwier. 1927.  _Manual of political economy Tr. by Ann S. Schwier_.  Macmillan, London, London and Basingstoke.
  * Stadler (1979) W. Stadler. 1979\.  A survey of multicriteria optimization or the vector maximum problem, part I: 1776-1960.  _Journal of Optimization Theory and Applications_ 29, 1 (1979), 1–52.  https://doi.org/10.1007/BF00932634
  * Wittneben (2009) Bettina B.F. Wittneben. 2009\.  Exxon is right: Let us re-examine our choice for a cap-and-trade system over a carbon tax.  _Energy Policy_ 37, 6 (2009), 2462–2464.  https://doi.org/10.1016/j.enpol.2009.01.029
  * Zhou et al. (2019) Dequn Zhou, Yunfei An, Donglan Zha, Fei Wu, and Qunwei Wang. 2019.  Would an increasing block carbon tax be better? A comparative study within the Stackelberg Game framework.  _Journal of Environmental Management_ 235, July 2018 (2019), 328–341.  https://doi.org/10.1016/j.jenvman.2019.01.082
  * Zitzler and Thiele (2006) Eckart Zitzler and Lothar Thiele. 2006.  Multiobjective optimization using evolutionary algorithms - A comparative case study.  (2006), 292–301.  https://doi.org/10.1007/bfb0056872