Estimation of Markovian-regime-switching models with independent regimes

TL;DR

This paper introduces new computational methods for estimating Markovian-regime-switching models with independent regimes, particularly applied to electricity prices, using augmented state techniques to improve efficiency.

Contribution

It develops forward, backward, and EM algorithms for independent regime MRS models by augmenting the hidden process with a counter for recent state visits.

Findings

Algorithms are computationally feasible for practical use.

Enhanced modeling of electricity prices with independent regimes.

Improved estimation accuracy over existing methods.

Abstract

Markovian-regime-switching (MRS) models are commonly used for modelling economic time series, including electricity prices where independent regime models are used, since they can more accurately and succinctly capture electricity price dynamics than dependent regime MRS models can. We can think of these independent regime MRS models for electricity prices as a collection of independent AR(1) processes, of which only one process is observed at each time; which is observed is determined by a (hidden) Markov chain. Here we develop novel, computationally feasible methods for MRS models with independent regimes including forward, backward and EM algorithms. The key idea is to augment the hidden process with a counter which records the time since the hidden Markov chain last visited each state that corresponding to an AR(1) process.

Tables

Table 1. Table 1: Distributions used to sample initial parameter values for the EM algorithm.

Parameter	${\alpha }_i$	${\varphi }_i$	${\sigma }_i^2$	$\mu$	$p_{11}$	$p_{22}$
Distribution	$U(-1,1)$	$U(-1,1)$	$U(0,4)$	$U(0,8)$	$U(0,1)$	$U(0,1)$

Table 2. Table 2: Proportion of terminating values which lie within 0.0001 of the optimal parameter value, as measured by the sup norm , | ⋅ | ∞ |\cdot|_{\infty} .

Simulation	1	2	3	4	5
Model 1	1	1	0.86	0.14	1
Model 2	0.96	1	1	1	1

Table 3. Table 3: BIC values of models M1-LN, M1-Gamma, M2-LN, and M2-Gamma, fitted to the SA electricity market data.

Model	M1-LN	M1-Gamma	M2-LN	M2-Gamma
BIC	15572	15681	15575	15706

Table 4. Table 4: Parameter estimates of model M1-LN, which has the lowest BIC, fitted to the SA electricity market data.

$B_t$	$\alpha$	$\varphi$	${\sigma }_1^2$		$p_{11}$
	-3.257	0.6830	213.26		0.9140
$S_t$	$q_3$	${\mu }_2$	${\sigma }_2^2$		$p_{22}$
	7.106	3.751	1.268		0.3945

Equations

$$X_t|\{R_t,X_{t-1},X_{t-2},\dots ,X_0\}\sim F^{R_t},$$

$$X_t|\{R_t=1,X_{t-1},X_{t-2},\dots ,X_0\}=0.6X_{t-1}+{\epsilon }_t^{\left(1\right)},$$

$$X_t|\{R_t=2,X_{t-1},X_{t-2},\dots ,X_0\}=1+0.9X_{t-1}+{\epsilon }_t^{\left(2\right)}.$$

$$B_t$$

$$S_t$$

$$X_t=\{\begin{array}{ll}{B}_t,& \text{ if }{R}_t=1,\\ S_t,& \text{ if }{R}_t=2.\end{array}$$

$$L\left(\mathbf{\theta }\right):=f_{\mathbf{X}}^{\mathbf{\theta }}\left(\mathbf{x}\right)=\sum _{\mathbf{R}\in {\mathcal{S}}^{T+1}}{f}_{\mathbf{X},\mathbf{R}}^{\mathbf{\theta }}\left(\mathbf{x},\mathbf{R}\right)=\sum _{\mathbf{R}\in {\mathcal{S}}^{T+1}}{f}_{\mathbf{X}|\mathbf{R}}^{\mathbf{\theta }}\left(\mathbf{x}|\mathbf{R}\right)f_{\mathbf{R}}^{\mathbf{\theta }}\left(\mathbf{R}\right),$$

$$f_{X_t|{\mathbf{X}}_{0:t-1}}^{\mathbf{\theta }}\left(x_t|{\mathbf{x}}_{0:t-1}\right)$$

$$=\sum _{i\in \mathcal{S}}{f}_{X_t|R_t,{\mathbf{X}}_{0:t-1}}^{\mathbf{\theta }}\left(x_t|i,{\mathbf{x}}_{0:t-1}\right)\sum _{j\in \mathcal{S}}{\mathbb{P}}^{\mathbf{\theta }}\left(R_{t-1}=j|{\mathbf{x}}_{0:t-1}\right)p_{ji},$$

$${\mathbb{P}}^{\mathbf{\theta }}\left(R_{t-1}=j|{\mathbf{x}}_{0:t-1}\right)$$

$$=\frac{f_{X_{t-1}|R_{t-1},{\mathbf{X}}_{0:t-2}}^{\mathbf{\theta }}\left(x_{t-1}|j,{\mathbf{x}}_{0:t-2}\right)\sum _{i\in \mathcal{S}}{p}_{ij}{\mathbb{P}}^{\mathbf{\theta }}\left(R_{t-2}=i|{\mathbf{x}}_{0:t-2}\right)}{f_{X_{t-1}|{\mathbf{X}}_{0:t-2}}^{\mathbf{\theta }}\left(x_{t-1}|{\mathbf{x}}_{0:t-2}\right)},$$

$$\log f_{\mathbf{X},\mathbf{R}}^{\mathbf{\theta }}\left(\mathbf{x},\mathbf{R}\right)$$

$$=\sum _{j\in \mathcal{S}}\mathbb{I}\left(R_0=j\right)\log f_{X_0|R_0}^{\mathbf{\theta }}\left(x_0|j\right)$$

$$+\sum _{t=1}^T\sum _{j\in \mathcal{S}}\mathbb{I}\left(R_t=j\right)\log f_{X_t|R_t,{\mathbf{X}}_{0:t-1}}^{\mathbf{\theta }}\left(x_t|j,{\mathbf{x}}_{0:t-1}\right)$$

$$+\sum _{i,j\in \mathcal{S}}{\eta }_{ij}\log p_{ij}+\sum _{j\in \mathcal{S}}\mathbb{I}\left(R_0=j\right)\log {\pi }_j,$$

$$Q\left(\mathbf{\theta },{\mathbf{\theta }}_n\right)$$

$$+\sum _{t=1}^T\sum _{j\in \mathcal{S}}{\mathbb{P}}^{{\mathbf{\theta }}_n}\left(R_t=j|{\mathbf{x}}_{0:T}\right)\log f_{X_t|R_t,{\mathbf{X}}_{0:t-1}}^{\mathbf{\theta }}\left(x_t|j,{\mathbf{x}}_{0:t-1}\right)$$

$$+\sum _{i,j\in \mathcal{S}}\mathbb{E}[{\eta }_{ij}|{\mathbf{x}}_{0:T};{\mathbf{\theta }}_n]\log p_{ij}+\sum _{j\in \mathcal{S}}{\mathbb{P}}^{{\mathbf{\theta }}_n}\left(R_0=j|{\mathbf{x}}_{0:T}\right)\log {\pi }_j,$$

$${\mathbb{P}}^{{\mathbf{\theta }}_n}\left(R_t=j,R_{t+1}=i|{\mathbf{x}}_{0:T}\right)$$

$${\mathbb{P}}^{{\mathbf{\theta }}_n}\left(R_t=j|{\mathbf{x}}_{0:T}\right)$$

$${\varphi }_j^{\left(n+1\right)}$$

$${\alpha }_j^{\left(n+1\right)}$$

$${\left({\sigma }_j^2\right)}^{\left(n+1\right)}$$

$$B_{1,t}^{\left(i\right)}$$

$$\text{mbox}and{B}_{2,t}^{\left(i\right)}$$

$$p_{ij}^{\left(n+1\right)}=\frac{\sum _{t=1}^T{\mathbb{P}}^{{\mathbf{\theta }}_n}\left(R_t=j|{\mathbf{x}}_{0:T}\right)\frac{p_{ij}^{\left(n\right)}{\mathbb{P}}^{{\mathbf{\theta }}_n}\left(R_{t-1}=i|{\mathbf{x}}_{0:t-1}\right)}{{\mathbb{P}}^{{\mathbf{\theta }}_n}\left(R_t=j|{\mathbf{x}}_{0:t-1}\right)}}{\sum _{t=1}^T{\mathbb{P}}^{{\mathbf{\theta }}_n}\left(R_{t-1}=i|{\mathbf{x}}_{0:T}\right)},$$

$$B_t^i:=\mathbb{I}\left(R_t=i\right)x_t+\mathbb{I}\left(R_t\ne i\right)\left({\alpha }_i+{\varphi }_i{B}_{t-1}^i+{\sigma }_i{\epsilon }_t^i\right).$$

$${\tilde{b}}_{t,i}^{\left(n\right)}$$

$$X_t=\{\begin{array}{ll}{B}_t,& \text{ if }{R}_t=1,\\ Y_t,& \text{ if }{R}_t=2,\end{array}$$

$$X_t=\{\begin{array}{ll}{B}_t^1,& \text{mbox}if{R}_t=1,\\ ⋮& \\ B_t^k,& \text{mbox}if{R}_t=k,\\ S_t^{k+1},& \text{mbox}if{R}_t=k+1,\\ ⋮& \\ S_t^M,& \text{mbox}if{R}_t=M,\end{array}$$

$$\{{\mathbf{H}}_t{\}}_{t\in \mathbb{N}}:=\{\left({\mathbf{N}}_t,R_t\right){\}}_{t\in \mathbb{N}}=\{\left(N_{t,1},\dots ,N_{t,k},R_t\right){\}}_{t\in \mathbb{N}},$$

$${\mathbb{P}}^{\mathbf{\theta }}\left({\mathbf{H}}_{t+1}=\left({\mathbf{n}}_{t+1},j\right)|{\mathbf{H}}_t=\left({\mathbf{n}}_t,i\right)\right)$$

$$=\{\begin{array}{ll}{p}_{ij}& \text{for }i\in {\mathcal{S}}_{AR}^c,j\in \mathcal{S},{\mathbf{n}}_{t+1}={\mathbf{n}}_t+1,\\ p_{ij}& \text{for }i\in {\mathcal{S}}_{AR},j\in \mathcal{S},{\mathbf{n}}_{t+1}={\mathbf{n}}_t^{\left(-i\right)}+1,\\ 0& \text{otherwise}.\end{array}$$

$$\mathbb{P}\left({\mathbf{H}}_0=\left(n_{0,1},\dots ,n_{0,k},j\right)\right)=\{\begin{array}{ll}{\pi }_j,& \text{for }j\in \mathcal{S},n_{0,i}=1,\text{ for all }i\in {\mathcal{S}}_{AR},\\ 0,& \text{otherwise}.\end{array}$$

$${\alpha }_{{\mathbf{n}}_t}^{\left(t\right)}\left(j\right):=f_{{\mathbf{H}}_t,{\mathbf{X}}_{0:t}}^{\mathbf{\theta }}\left(\left({\mathbf{n}}_t,j\right),{\mathbf{x}}_{0:t}\right)$$

$${\alpha }_{{\mathbf{n}}_0}^{\left(0\right)}\left(j\right)=f_{X_0|{\mathbf{H}}_0}^{\mathbf{\theta }}\left(x_0|\left({\mathbf{n}}_0,j\right)\right){\mathbb{P}}^{\mathbf{\theta }}\left({\mathbf{H}}_0=\left({\mathbf{n}}_0,j\right)\right).$$

$${\alpha }_{{\mathbf{n}}_t}^{\left(t\right)}\left(j\right)$$

$$L\left(\mathbf{\theta }\right)=\sum _{j\in \mathcal{S}}\sum _{{\mathbf{n}}_T\in {\mathcal{S}}^{\left(T\right)}}{\alpha }_{{\mathbf{n}}_T}^{\left(T\right)}\left(j\right).$$

Full text

Keywords: Electricity price model, forward-backward algorithm, hidden Markov model, Markov-switching time series

1 Introduction

A commonly used model for economic time series is the Markovian-regime-switching (MRS) model whereby multiple stochastic processes are interweaved by a Markov chain. The general idea is that there exist multiple regimes underlying the observation process, and depending on which regime the system is in, different characteristics are displayed. For example, for stock prices we could suppose that there is a bull regime where prices trend upward and are comparatively non-volatile, and a bear regime where prices trend downward and are relatively volatile. Our motivating application is electricity prices where it is common to model prices with an MRS model. Due to the fact that electricity cannot currently be stored efficiently, electricity prices show characteristics not seen in typical commodity markets, for example mean reversion, prices spikes, drops and negative prices. MRS models are able to capture these behaviours, and have been popular tools for modelling randomness in electricity markets. Typically MRS models with two or three regimes are used [9, 8, 13, 14]: a base regime where prices are relatively non-volatile, a spike regime where prices are volatile and high, and sometimes a drop regime is included, where prices are volatile and low. More broadly, models with Markovian switching find application in biology [1], weather modelling [23, 27], speech recognition [21] and more.

The earliest applications of MRS models for electricity prices are [9] and [8]. Their models specify that prices decay back to base levels following a spike according to an autoregressive process of order 1 (AR(1)). However, this is not consistent with observations from the market where a more immediate return to base levels is observed [16, 2]. For this reason, [14] introduce a three-regime MRS model which separates the behaviour of base and spike prices. They use one regime to capture base prices, one regime to capture spikes, and one regime to return prices to base levels following a spike. Motivated by the need to capture the distinct and abrupt price spikes in electricity markets, it has become popular to specify MRS models with independent regimes, which are able to capture this behaviour without the addition of the extra regime required by [14]. We may think of a dependent regime MRS models as a single process, where the dynamics of the process is governed by the hidden Markov chain, whereas we may think of an independent regime MRS model as a collection of independent AR(1) processes, and at each time $t$, the hidden Markov chain chooses which process is observed.

We say a model has independent regimes if, given the hidden regime sequence, the observations generated from each regime are independent of observations generated from any other regime; and we say a model has dependent regimes otherwise. Independent regime MRS models were introduced in [13] and have since been popular [25, 17, 24]. In these models, typically at least one regime is specified as an AR(1) process. AR(1) processes have a dependence structure between values at successive times which, coupled with the independent regimes assumption in an MRS model, complicates the dependence structure between an observation at time $t$ and all prior observations – the dependence between prices is governed by the hidden regime process and is therefore random. It is for this reason that the forward, backward and expectation-maximisation (EM) algorithms for traditional (dependent regime) MRS models or hidden Markov models do not apply. For each time $t$, the forward algorithm evaluates the probabilities that the hidden process is in each regime given the observed values up to time $t$. The backward algorithm then uses the output of the forward algorithm to calculate the probabilities that the hidden process is in each regime given all observations. The EM algorithm is an iterative optimisation algorithm, iterating between an E-step and an M-step, used to find the maximum likelihood estimates of parameters for models with missing data or latent variables – the E-step is computed by the backward algorithm.

For simplicity, in this work, we focus on independent regime MRS models with AR(1) and i.i.d. regimes only, since these are the type of models used in the electricity price modelling literature. However, we believe the methods developed here are more general and apply to Markov-switching processes with regimes that are discrete-time Markov chains generally, and can also be extended to more general autoregressive processes. The most popular method of inference for the models used for electricity pricing is an approximation to the EM algorithm introduced by [17], which we show can be unreliable (see [22] also). Here, a novel, computationally feasible, and exact likelihood-based framework to solve this problem is developed. This work is related to the forward, backward, and EM algorithms for traditional MRS models, and, more closely, to the same algorithms for hidden semi-Markov models, where the idea of augmenting the hidden process with a counter is also used [26]. The novel algorithms have complexity $\mathcal{O}\left(M^{2}T^{k+1}k^{k}\right)$ where $M<\infty$ is the total number of regimes in the model, $T$ is the length of the observed data set, and $k$ is the number of AR(1) regimes in the model. As $\mathcal{O}\left(M^{2}T^{k+1}k^{k}\right)$ may be impractically large depending on the values of $T$ and $k$, we also present an approximation to our algorithm that is $\mathcal{O}\left(M^{2}TD^{k}k^{k}\right)$ where $D>0$ is the maximum length of the memory of the AR(1) processes, that is, we specify that the observations $x_{t}$ may depend only on $x_{t-\ell}$ for $\ell\in\{1,2,...,D\}$, and is independent of any $x_{t-\ell}$ for $\ell>D$.

This paper is structured as follows. We formally introduce the MRS model, in particular the independent-regime type in Section 1.1, and discuss the approximate parameter inference algorithm of [17] in Section 1.2. A novel forward algorithm, which is used to evaluate the likelihood and filtered state probabilities for these models, is presented in Section 2. Using the outputs of the forward algorithm we develop a novel backward algorithm in Section 3. The backward algorithm is used to evaluate the smoothed state probabilities, which are applied in Section 4 to construct an EM algorithm. We introduce the truncated approximations of our algorithms at the end of Section 4. Section 5 provides simulation evidence that our EM algorithm is consistent and that the truncation approximations are reasonable, while in Section 6 we apply our algorithms to estimate MRS models for the South Australian wholesale electricity market. Finally, we make concluding remarks in Section 7.

1.1 MRS models: A brief introduction

An MRS model is built from two pieces, an unobservable regime sequence, $\{R_{t}\}_{t\in\mathbb{N}}$, which is a finite-state Markov Chain, and an observation sequence, $\{X_{t}\}_{t\in\mathbb{N}}$. Let us denote the state space of the hidden regime sequence as $\mathcal{S}=\{1,2,\dots,M<\infty\}$, and the transition matrix as $P=[p_{ij}]_{i,j\in\mathcal{S}}$.

The simplest MRS model is the hidden Markov model (HMM) where observations $X_{t}$ take values in a discrete set, and $X_{t}$ is independent of $X_{t-1},\dots,X_{0}$ and $X_{t+1},X_{t+2},\dots$ given the regime at time $t$, $R_{t}$. In general, MRS models are specified in terms of distributions that allow dependence on past observations, given the current regime. That is, the model defines distributions,

[TABLE]

for some distribution $F^{R_{t}}$. The MRS model, as introduced by [10, 11], specifies that $X_{t}|\{R_{t},X_{t-1},X_{t-2},\dots,X_{0}\}$ follows some time-series model (an autoregressive process of order $p$, for example) with dependence on a finite number of past observations, but not on $R_{0},..,R_{t-1}$. That is, the dependence structure does not take into account which regime the past observations belong to. For example, the following is a dependent regime MRS model.

Example 1** (An MRS model with dependent regimes).**

Let $\mathcal{S}=\{1,2\}$, and $p_{11}=p_{22}=0.9$, and specify

[TABLE]

for $\varepsilon_{t}^{\left(i\right)}\sim$ i.i.d. N(0,1) for $i=1,2$. So $X_{t}$ follows AR(1) dynamics in both Regimes 1 and 2. This is a dependent-regime MRS model since $X_{t}$ depends on $X_{t-1}$ regardless of which regime the lagged observation, $X_{t-1}$, came from. Figure 1 (Left) shows a simulation of this model.

In this paper we relax the assumption that the current observation $X_{t}$ is conditionally independent of $R_{0},..,R_{t-1}$, in order to increase flexibility in these models. In particular, we consider models where, given $R_{t}=i$, $X_{t}$ depends only on lagged values from Regime $i$, thus the dependence structure is random as it is a function of $R_{0},\dots,R_{t-1}$. The following is an example of an independent regime MRS model.

Example 2** (An MRS model with independent regimes).**

Let $\mathcal{S}=\{1,2\}$, and $p_{11}=p_{22}=0.9$, and define the following AR(1) processes

[TABLE]

where $\varepsilon_{t}^{B}$ and $\varepsilon_{t}^{S}$ are sequences of i.i.d. N(0,1) random variables. Then, construct the MRS model as follows

[TABLE]

Figure 1 (Right) shows a simulation of this model.

A precise definition of dependent and independent regime models is the following. Define the sets $\mathcal{A}_{i}:=\{t\in\mathbb{N}:R_{t}=i\}$, $i\in\mathcal{S}$. We say that a model has independent regimes if, given the regime sequence, the sets $\{X_{t}:t\in\mathcal{A}_{i}\}$, $i\in\mathcal{S}$, are independent. Otherwise, it is a dependent regime model.

1.2 Existing methods

For the simplest form of MRS model, the HMM, likelihood evaluation and maximisation algorithms were first developed in a series of papers, [5], [4], and [6], and subsequent work on MRS models is typically closely related to this. The first algorithms for the more general dependent regime MRS model were presented by [10, 11]. The main issue for maximum likelihood estimation of models with hidden regimes is that the regime sequence is unobserved, thus to naively evaluate the likelihood requires calculation of the marginal distribution

[TABLE]

where $f_{\boldsymbol{Y}}^{\boldsymbol{\theta}}\left(\cdot\right)$ denotes the distribution function of a random vector $\boldsymbol{Y}$ with parameters $\boldsymbol{\theta}$, $\boldsymbol{x}=\left(x_{0},\dots,x_{T}\right)$ is a sequence of observed values, and $\mathcal{S}^{T+1}$ is the space of all possible regime sequences of length $T+1$, $\boldsymbol{R}=\left(R_{0},...,R_{T}\right)$. The number of sequences in $\mathcal{S}^{T+1}$ is $M^{T+1}$ which, for most realistic datasets, is computationally infeasible to evaluate in this form. In the context of HMMs, the sum (1) is made computationally feasible by the forward algorithm [6], and the maximisation of the likelihood is commonly performed via the Baum-Welch algorithm [6], which is a specific case of the EM algorithm [7] and uses the backward algorithm [6].

The works of [10, 11] extend the methods for HMMs to MRS models with dependent regimes by adapting the forward algorithm, developing a new algorithm to replace the backward algorithm and constructing an EM algorithm. [20] refines the work of Hamilton, developing a more efficient implementation of Hamilton’s smoothing algorithm. Kim’s algorithm is similar to the backward algorithm for HMMs.

Relevant to this paper, [17] extend Hamilton’s work and develop an approximate algorithm for MRS models with independent regimes, which we label the EM-like algorithm since it resembles Hamilton’s EM algorithm. However, it is not an example of the EM algorithm and so none of the EM theory holds. In Sections 1.2.1–1.2.3, we briefly review the work of [11], [20] and [17] in order to provide the motivation and background for our work.

1.2.1 Likelihood evaluation for dependent-regime models: The forward algorithm

Define $\boldsymbol{x}_{r:s}=\left(x_{r},x_{r+1},\dots,x_{s}\right)$ for $r\leq s$ and write the likelihood as $L\left(\boldsymbol{\theta}\right)=f^{\boldsymbol{\theta}}_{X_{0}}\left(x_{0}\right)\prod\limits_{t=1}^{T}f^{\boldsymbol{\theta}}_{X_{t}|\boldsymbol{X}_{0:t-1}}\left(x_{t}|\boldsymbol{x}_{0:t-1}\right).$ The forward algorithm [11] calculates $f^{\boldsymbol{\theta}}_{X_{0}}$ and $f^{\boldsymbol{\theta}}_{X_{t}|\boldsymbol{X}_{0:t-1}}$ for $t=1,2,\dots,T,$ from which it is straightforward to calculate the likelihood or loglikelihood:

Algorithm 1: The forward algorithm [11]

1. Step 1.

Initialise the algorithm with values $\mathbb{P}^{\boldsymbol{\theta}}\left(R_{0}=i\right)=\pi_{i}$, which may be assumed to be known a priori, or left as parameters to be inferred. 2. Step 2.

The term, $f^{\boldsymbol{\theta}}_{X_{0}}\left(x_{0}\right)$, is calculated as $f^{\boldsymbol{\theta}}_{X_{0}}\left(x_{0}\right)=\sum\limits_{i\in\mathcal{S}}f^{\boldsymbol{\theta}}_{X_{0}|R_{0}}\left(x_{0}|i\right)\mathbb{P}^{\boldsymbol{\theta}}\left(R_{0}=i\right),$ where the density $f^{\boldsymbol{\theta}}_{X_{0}|R_{0}}$ is known from the model specification. 3. Step 3.

For $t=1,\dots,T$:

[TABLE]

where $f^{\boldsymbol{\theta}}_{X_{t}|R_{t},\boldsymbol{X}_{0:t-1}}$ is also known from the model specification. Furthermore, the probabilities $\mathbb{P}^{\boldsymbol{\theta}}\left(R_{t-1}=j|\boldsymbol{x}_{0:t-1}\right)$ can be calculated using Bayes’ Theorem,

[TABLE]

for $t\geq 2$. These are known as the forward/filtered probabilities.

The quantities $\mathbb{P}^{\boldsymbol{\theta}}\left(R_{t}=i|\boldsymbol{x}_{0:t-1}\right)=\sum\limits_{j\in\mathcal{S}}p_{ji}\mathbb{P}^{\boldsymbol{\theta}}\left(R_{t-1}=j|\boldsymbol{x}_{0:t-1}\right)$ are known as the prediction probabilities. In some applications, the forward and prediction probabilities may be quantities of interest in their own right, and they also appear as inputs to the backward algorithm (see Section 1.2.2).

1.2.2 Maximum likelihood for dependent-regime models: The EM algorithm

Hamilton’s forward algorithm [11] is a computationally feasible way to evaluate the loglikelihood, from which it is possible to use black-box optimisation methods to find the MLEs. However, it is common to use the EM algorithm [7] instead, particularly when the E-step and M-step of the algorithm are available in closed form. The EM algorithm proceeds by iterating between the E-step, constructing the function $Q\left(\boldsymbol{\theta},\boldsymbol{\theta}_{n}\right)=\mathbb{E}\left[\log f^{\boldsymbol{\theta}}_{\boldsymbol{X},\boldsymbol{R}}\left(\boldsymbol{x},\boldsymbol{R}\right)\mid\boldsymbol{x};\boldsymbol{\theta}_{n}\right]$, and the M-step, maximising $Q$ with respect to $\boldsymbol{\theta}\in\Theta$, where $\Theta$ is the parameter space. This results in a sequence $\{\boldsymbol{\theta}_{n}\}_{n\in\mathbb{N}}$ that converges to a local maximiser of the loglikelihood.

The EM algorithm for MRS models with dependent regimes proceeds as follows [11]. Define the random variable $\eta_{ij}$ as the number of transitions from state $i$ to state $j$ in the sequence $\boldsymbol{R}=\left(R_{0},\dots,R_{T}\right)$ and let $\mathbb{I}\left(\cdot\right)$ be the indicator function. The joint log-density of $\boldsymbol{x}$ and $\boldsymbol{R}$ can be written as

[TABLE]

where $\pi_{j}$ denotes $\mathbb{P}^{\theta}\left(R_{0}=i\right)$. In the $n$th iteration, $n\geq 0$, for the E-step, taking the conditional expectation given parameters $\boldsymbol{\theta}_{n}$ and observed values $\boldsymbol{x}_{0:T}$ yields

[TABLE]

where the expectation $\mathbb{E}[\eta_{ij}|\boldsymbol{x}_{0:T};\boldsymbol{\theta}_{n}]=\sum\limits_{t=1}^{T}\mathbb{P}^{\boldsymbol{\theta}_{n}}\left(R_{t}=j,R_{t-1}=i|\boldsymbol{x}_{0:T}\right)$. The densities $f^{\boldsymbol{\theta}}_{X_{0}|R_{0}}$ and $f^{\boldsymbol{\theta}}_{X_{t}|R_{t},\boldsymbol{X}_{0:t-1}}$, are given by the model specification.

The smoothed probabilities, $\mathbb{P}^{\boldsymbol{\theta}_{n}}\left(R_{t}=j|\boldsymbol{x}_{0:T}\right)$ and $\mathbb{P}^{\boldsymbol{\theta}_{n}}\left(R_{t}=j,R_{t-1}=i\mid\boldsymbol{x}_{0:T}\right)$, required to construct $Q$ are obtained using a backward recursion after running the forward algorithm with parameters $\boldsymbol{\theta}_{n}$, and storing the forward and prediction probabilities. Developed by [20], this backward recursion is in Algorithm 2, below.

Algorithm 2: The backward algorithm [20]

1. Step 1.

Evaluate $\mathbb{P}^{\boldsymbol{\theta}_{n}}\left(R_{T}=j|\boldsymbol{x}_{0:T}\right)$ using the forward algorithm (Algorithm 1). 2. Step 2.

For $t=T-1,\dots,0,$

[TABLE]

where $p_{ij}^{\left(n\right)}$ means the $p_{ij}$ parameter under $\boldsymbol{\theta}_{n}$.

After executing Kim’s backward algorithm, we can construct the function $Q$. In the M-step, the maximisers of $Q\left(\cdot,\boldsymbol{\theta}_{n}\right)$ are found. In the dependent regime model, if the process is in Regime $j$ at time $t$, then the observations evolve according to $X_{t}=\alpha_{j}+\phi_{j}X_{t-1}+\sigma_{j}\varepsilon_{t}$ where $\alpha_{j},$ $\phi_{j}$ and $\sigma_{j}$ are parameters, and $\{\varepsilon_{t}\}\sim\mbox{N}\left(0,1\right)$. Recall that, for the dependent regime model, $X_{t-1}$ is the last observed value (regardless of which regime generated it). The maximiser of $Q$ at the $\left(n+1\right)$th iteration of the EM algorithm, $\boldsymbol{\theta}^{n+1}$, for $n\geq 0$, is the following system of equations [11, 17]:

[TABLE]

In general, the switching probabilities are updated using the following [20]

[TABLE]

which rely on the smoothed, forward and prediction probabilities. For i.i.d. regimes the M-step can often be derived analytically; however, such expressions are not required for our discussion since there is no dependence on lagged values in these regimes and so they are omitted.

Thus, we implement the EM algorithm by initialising it with a guess of the true parameters, then alternating between the forward and backward algorithms (the E-step) and calculating the maximisers of $Q$ (the M-step). The algorithm terminates when the step size is below a prespecified tolerance, i.e. $|\boldsymbol{\theta}_{n+1}-\boldsymbol{\theta}_{n}|_{\infty}<e$ where $e$ is some small tolerance.

1.2.3 Approximate maximum likelihood for independent-regime models: The EM-like algorithm

For independent-regime MRS models the EM algorithm is computationally infeasible if the densities $f^{\boldsymbol{\theta}}_{X_{t}|R_{t},\boldsymbol{X}_{0:t-1}}$, $R_{t}\in\mathcal{S}$, $t=1,\dots,T$, are computed naively, as this is a $\mathcal{O}\left(M^{t}\right)$ calculation for each $t=1,...,T$, where $M<\infty$ is the number of regimes in the model. On the other hand, if we use our proposed foward algorithm, introduced in Section 2), calculation of these densities is $\mathcal{O}\left(t^{k+1}\right)$ for each $t=1,...,T$, which may be computationally feasible when $T$ and $k$ are not too large.

Developed by [17], the EM-like algorithm is an approximation to the EM algorithm. For independent regime MRS models, the EM-like algorithm overcomes the problem of computational infeasibility by replacing lagged values for Regime $i$ (assuming this is an AR(1) regime) with approximations, $\tilde{b}_{t-1,i}^{\left(n\right)}$. These are described as the expectations [17], $\mathbb{E}[B^{i}_{t}|\boldsymbol{x}_{0:t};\boldsymbol{\theta}_{n}],$ where

[TABLE]

Simply put, wherever $x_{t-1}$ appears verbatim in an expression related to Regime $i$ in the EM algorithm in Sections 1.2.1-1.2.2, it is replaced with $\tilde{b}_{t-1,i}^{\left(n\right)}$ at the $n^{\text{th}}$ iteration. The $\tilde{b}_{t,i}^{\left(n\right)}$ are calculated recursively as

[TABLE]

where $\tilde{\mathbb{P}}^{\boldsymbol{\theta}_{n}}\left(R_{t}=i|\boldsymbol{x}_{0:t}\right)$ and $\tilde{\mathbb{P}}^{\boldsymbol{\theta}_{n}}\left(R_{t}\neq i|\boldsymbol{x}_{0:t}\right)$ are given by the forward algorithm which is part of the EM-like procedure. Janczura and Weron, [17], conduct simulation studies and show that this algorithm seems to work well for the datasets they generate. However, no theoretical results are available that show convergence of, or error bounds for, the EM-like algorithm; in particular, there is no guarantee that the parameter estimates produced by the EM-like algorithm are consistent. In contrasts, our algorithms rest on the theory of the EM algorithm.

We can construct examples of independent-regime MRS models where the EM-like algorithm fails to get close to the true parameter values.

Example 3**.**

Consider the following independent-regime MRS model,

[TABLE]

where $B_{t}$ is an AR(1) process, $B_{t}=0.95B_{t-1}+\sqrt{0.2}\varepsilon_{t},$ with $\{\varepsilon_{t}\}$ being a sequence of i.i.d. N$\left(0,1\right)$ random variables, $Y_{t}$ is an i.i.d. sequence of N$\left(2,1\right)$ random variables, and $\{R_{t}\}_{t\in\mathbb{N}}$ is a Markov chain with state space $\mathcal{S}=\{1,2\}$, transition matrix entries $p_{11}=0.5$ and $p_{22}=0.8$, and initial probability distribution $\left(1,0\right)$, so the process always starts in Regime 1. We simulated 20 realisations of length $T=2000$ from this model and used the EM-like algorithm to try to recover the true parameters. To give the algorithm the best chance of converging to the true parameters, we initialise the EM-like algorithm at the true parameter values. The parameters recovered by the EM-like algorithm are summarised in Figure 2.

For comparison, the MLEs obtained using our EM algorithm (Section 4) are also shown. Notice in Figure 2 that the EM-like algorithm performs poorly, while our exact method performs much better.

2 A novel forward algorithm

The general idea of the novel algorithms presented in this paper is to augment the hidden Markov chain with counters that record the last time each AR(1) regime was visited. This augmented process is a Markov chain, and similar arguments to those used to construct the forward-backward algorithm for MRS models with dependent regimes can be used to construct a forward and backward algorithms for these models. Our methods are related to the forward and backward algorithms for hidden semi-Markov models, where the hidden process is also augmented with a counter and the augmented hidden process is a Markov chain [26]. Though similar, the algorithms for hidden semi-Markov models (HSMMs) are not applicable to the models considered here since for HSMMs the counter counts the number of time steps since the last change of regime, whereas here the counters count the number of transitions since the last visit to a regime. Furthermore, in HSMMs, the counters do not appear in the conditional densities of the observations, as they do here.

2.1 The augmented hidden Markov chain

For the following, suppose that the first $k$ states, $\mathcal{S}_{AR}=\{1,\dotsc,k<M\}$, correspond to AR(1) processes and all other regimes are i.i.d. That is, we have the following independent regime MRS model

[TABLE]

where $B_{t}^{i}=\alpha_{i}+\phi_{i}B^{i}_{t-1}+\sigma_{i}\varepsilon_{t}^{i}$ are AR(1) and $S_{t}^{j}$ are i.i.d. Our arguments also hold for $k=M$ with only slight modification, but here we treat the case $k<M$ only, since these are the types of models relevant to our application.

Now, define another Markov chain

[TABLE]

where $N_{t,j}\in\mathbb{N}_{+}$ counts the number of time steps since the process $\{R_{t}\}$ was last in Regime $j$ before time $t$, for each AR(1) regime $j=1,\dots,k$. When there is no time $\tau\in\{0,1,\dots,t-1\}$ with $R_{\tau}=j$, then we set $N_{t,j}=t+1$. Thus the augmented Markov chain $\{\boldsymbol{H}_{t}\}_{t\in\mathbb{N}}$ lives on the state space $\mathbb{N}_{+}^{k}\times\mathcal{S}.$

To describe the transitions of the Markov chain $\{\boldsymbol{H}_{t}\}$, let $\boldsymbol{n}:=\left(n_{1},\dots,n_{k}\right)\in~{}\mathbb{N}_{+}^{k}$ be an arbitrary vector of counters, with $n_{r}\neq n_{s}$ for $r\neq s$ or, at time $t$, $n_{r}=n_{s}=t+1$ (it is possible that neither state $r$ nor state $s$ have been visited by time $t$). Also, define $\boldsymbol{1}$ to be a row vector of ones of length $k$, $\boldsymbol{e}_{i}$ to be a row vector of length $k$ with all entries being 0 except the $i^{\text{th}}$ entry which is 1, and $\boldsymbol{n}^{\left(-i\right)}:=\boldsymbol{n}-n_{i}\boldsymbol{e}_{i}=\left(n_{1},\dots,n_{i-1},0,n_{i+1},\dots,n_{k}\right).$

The transition probabilities of $\{\boldsymbol{H}_{t}\}$ are

[TABLE]

In words, when the current state is $\boldsymbol{H}_{t}=\left(\boldsymbol{n}_{t},i\right)$ and $i$ is not an AR(1) regime (so there is no counter associated with state $i$), then, at time ${t+1}$, $R_{t}$ transitions to state $j$ with probability $p_{ij}$ and all the counters are advanced by $1$ to $\boldsymbol{n}_{t+1}=\boldsymbol{n}_{t}+\boldsymbol{1}$, since there has been one more time step since $\{R_{t}\}$ was last in any state with a counter (any state in $\mathcal{S}_{AR}$). When the current state is $\boldsymbol{H}_{t}=\left(\boldsymbol{n}_{t},i\right)$, where $i$ is an AR(1) regime, then $R_{t}$ transitions to any state $j\in\mathcal{S}$ with probability $p_{ij}$, the counter for Regime $i$, $n_{{t+1},i}$, is set to 1, since the last time in state $i$ was $t$, and all other counters are advanced by $1$. All other transition probabilities for $\{\boldsymbol{H}_{t}\}$ are 0.

The state space of $\{\boldsymbol{H}_{t}\}$ is countably infinite. However, due to the way $\{\boldsymbol{H}_{t}\}_{t\in\mathbb{N}}$ is initialised and evolves, many states are inaccessible for $\{\boldsymbol{H}_{t}\}$, $t\in\{0,1,\dots,T<\infty\}$, and this makes our algorithm computationally feasible. Specifically, we suppose that the Markov chain $\{\boldsymbol{H}_{t}\}$ is initialised with the probability distribution

[TABLE]

The distribution $\boldsymbol{\pi}:=\left(\pi_{1},\dots,\pi_{M}\right)$ can be any proper probability distribution. However, in line with existing algorithms for dependent regime MRS models, it can be either the stationary distribution of $\{R_{t}\}$, or a point mass on a single state, or, when used as part of the EM algorithm, the probabilities $\mathbb{P}^{\boldsymbol{\theta}_{n}}\left(\boldsymbol{H}_{0}|\boldsymbol{x}_{0:T}\right)$ calculated at the previous iteration of the EM algorithm. The following lemma gives all the states that $\{\boldsymbol{H}_{t}\}$ can be in at time $t>0$.

Lemma 4**.**

Define $\mathcal{S}^{\left(0\right)}:=\{\boldsymbol{1}\}$ as a vector of $1$’s of length $k$ where $k$ is the number of AR(1) regimes. For each $t=1,2,\dots,T$, let $\mathcal{S}^{\left(t\right)}$ be the set of all vectors $\boldsymbol{n}_{t}:=\left(n_{t,1},\dots,n_{t,k}\right)$ such that, for $j,m\in\mathcal{S}_{AR}$,

1. (i)

$n_{t,j}\in\{1,2,\dots,t+1\}$, 2. (ii)

there are at most $\min\left(t,k\right)$ elements of $\boldsymbol{n}_{t}$ with $n_{t,j}\neq t+1$, 3. (iii)

$n_{t,j}\neq n_{t,m}$* for all $j\neq m$, unless $n_{t,j}=n_{t,m}=t+1$.*

Given $\{\boldsymbol{H}_{t}\}$ is initialised with the distribution in Equation (7), it is possible for $\boldsymbol{H}_{t}$ to reach states $\left(\boldsymbol{n}_{t},i\right)$ where $\boldsymbol{n}_{t}\in\mathcal{S}^{\left(t\right)}$ and $i\in\mathcal{S}$ only. The cardinality of $\mathcal{S}^{\left(t\right)}$ is $|\mathcal{S}^{\left(t\right)}|=\sum\limits_{m=0}^{\min\left(t,k\right)}\binom{t}{m}\binom{k}{m}m!.$

Proof.

First, we explain why $\mathcal{S}^{\left(t\right)}$ contains all possible values of the counters of $\boldsymbol{H}_{t}$ for $t\geq 0$.

At time $t=0$ the chain, $\{\boldsymbol{H}_{t}\}=\{\left(\boldsymbol{N}_{t},R_{t}\right)\}$, is initialised with the distribution in Equation (7), so $\mathcal{S}^{\left(0\right)}:=\{\boldsymbol{1}\}.$

At time $t>0$, either $\{R_{t}\}$ has never visited state $j\in\mathcal{S}_{AR}$, in which case $n_{t,j}=t+1$, or $\{R_{t}\}$ last visited $j$ at time $t_{j}$, in which case $n_{t,j}=t-t_{j}\in\{1,2,\dots,t\}$; this is part (i). Since the process $\{R_{t}\}$ can only be in one regime at a time, it follows that $n_{t,j}\neq n_{t,m}$ for $j\neq m$ (unless $n_{t,j}=n_{t,m}=t+1$), which is part (iii) of the definition. Also, at time $t$, the regime chain $\{R_{t}\}$ could only possibly have visited $\min\left(t,k\right)$ possible states; this is part (ii) of the definition.

Now, to prove the cardinality of $\mathcal{S}^{\left(t\right)}$. The elements of $\mathcal{S}^{\left(t\right)}$ are of the form $\left(n_{1},\dots,n_{k}\right)$. At time $t$, let $m$ be the possible number of counters that are not equal to $t+1$, so $m$ is an element of $\{0,1,\dots,\min\left(t,k\right)\}$. For each$m\in\{0,1,\dots,\min\left(t,k\right)\}$, there are $\binom{k}{m}$ ways of choosing which $m$ of the $k$ counters are not equal to $t+1$. Next, each counter takes a distinct value in $\{1,\dots,t\}$, so there are $\binom{t}{m}$ ways of choosing the value of the $m$ counters. There are $m!$ possible permutations to allocate the chosen values to the counters. So, in total there are $\sum\limits_{m=0}^{\min\left(t,k\right)}\binom{t}{m}\binom{k}{m}m!$ elements in $\mathcal{S}^{\left(t\right)}$. ∎

Lemma 4 says that if $\boldsymbol{n}_{t}\notin\mathcal{S}^{\left(t\right)}$ then $\mathbb{P}^{\boldsymbol{\theta}}\left(\boldsymbol{H}_{t}=\left(\boldsymbol{n}_{t},j\right)\right)=0$ for any $j\in\mathcal{S}$. Therefore the elements of the set $\mathcal{S}^{\left(t\right)}$ partition the space of all counters that the process $\{\boldsymbol{H}_{t}\}$ has positive probability of reaching. Thus, for any (measurable) set $A$ and any $t$, the law of total probability can be applied as $\mathbb{P}^{\boldsymbol{\theta}}\left(A\right)=\sum\limits_{\boldsymbol{n}_{t}\in\mathcal{S}^{\left(t\right)}}\sum\limits_{j\in\mathcal{S}}\mathbb{P}^{\boldsymbol{\theta}}\left(\boldsymbol{H}_{t}=\left(\boldsymbol{n}_{t},j\right),A\right).$ We will use this fact multiple times to construct the forward algorithm.

2.2 Constructing the forward algorithm

The forward algorithm is multi-purpose. It can be used to evaluate the likelihood, and also to evaluate the filtered and prediction probabilities, which in turn are inputs to the backward algorithm. For clarity of exposition, we first present a simple, but impractical due to underflow, algorithm (Lemma 5) to calculate the likelihood for independent regime MRS models, then address the underflow issue later with a normalised version of the algorithm (Lemma 7). Define

[TABLE]

for $t=0,1,\dots,T$, $\boldsymbol{n}_{t}\in\mathcal{S}^{\left(t\right)}$, and $j\in\mathcal{S}$.

Lemma 5** (A simple forward algorithm).**

First, for $j\in\mathcal{S}$ calculate

[TABLE]

Then for $t=1,2,\dots,T$, $\boldsymbol{n}_{t}\in\mathcal{S}^{\left(t\right)}$, $j\in\mathcal{S}$, calculate

[TABLE]

Then the likelihood is given by

[TABLE]

Proof.

First, from the law of total probability, we have

[TABLE]

Now, if any $n_{t,\ell}=1$ for $\ell\in\{1,\dots,k\}$, then it must be that $R_{t-1}=\ell$ and $\boldsymbol{n}_{t-1}=\boldsymbol{n}_{t}-\boldsymbol{1}+m\boldsymbol{e}_{\ell}$ for some $m\in\{1,\dots,t\}$. Thus, in this case, the double sum in Equation (11) simplifies to $\sum\limits_{m=1}^{t}f^{\boldsymbol{\theta}}_{\boldsymbol{H}_{t-1},\boldsymbol{H}_{t},\boldsymbol{X}_{0:t}}\left(\left(\boldsymbol{n}_{t}-\boldsymbol{1}+m\boldsymbol{e}_{\ell},\ell\right),\left(\boldsymbol{n}_{t},j\right),\boldsymbol{x}_{0:t}\right).$ Otherwise, all elements of $\boldsymbol{n}_{t}$ are greater than 1, in which case $R_{t-1}\notin\mathcal{S}_{AR}$ and $\boldsymbol{n}_{t-1}=\boldsymbol{n}_{t}-\boldsymbol{1}$, so the double sum in Equation (11) simplifies to

[TABLE]

For both cases the following arguments are the same, so for notational convenience we will use $\boldsymbol{n}$ to be either $\boldsymbol{n}_{t}-\boldsymbol{1}$ when $n_{t,\ell}>1$ for all $\ell\in\mathcal{S}_{AR}$, or $\boldsymbol{n}_{t}-\boldsymbol{1}+m\boldsymbol{e}_{\ell}$ when $n_{t,\ell}=1$ for some $\ell\in\mathcal{S}_{AR}$.

We can write the summands as

[TABLE]

and Equation (9) follows after noting that

[TABLE]

and from the definition of $\alpha_{\boldsymbol{n}}^{\left(t-1\right)}\left(i\right)$. Equation (10) is just an application of the law of total probability. ∎

Lemma 6**.**

The complexity of the simple forward algorithm, as given by the total number of multiplications, is less than $\mathcal{O}\left(M^{2}T^{k+1}k^{k}\right)$.

Proof.

For $t=0$, calculating $\alpha^{\left(0\right)}_{\boldsymbol{n}_{0}}\left(j\right)$ by (8) for all $j\in\mathcal{S}$ requires $M$ multiplications in total. For each $t\in\{1,\dots,T\}$, first consider the case where $n_{t,\ell}\neq 1$ for all $\ell\in\mathcal{S}_{AR}$. Fix $t$. Noting that

[TABLE]

then there are $|\mathcal{S}^{\left(t-1\right)}|\left(M-k\right)M$ multiplications to calculate all the necessary $p_{ij}\alpha_{\boldsymbol{n}_{t}-\boldsymbol{1}}^{\left(t-1\right)}\left(i\right)$ terms since $i\in\mathcal{S}_{AR}^{c}$, $j\in\mathcal{S}$ and $\boldsymbol{n}_{t}-\boldsymbol{1}\in\mathcal{S}^{\left(t-1\right)}$. The sum then collapses this to $|\mathcal{S}^{\left(t-1\right)}|M$ terms, each of which is then multiplied by $f^{\boldsymbol{\theta}}_{X_{t}|\boldsymbol{H}_{t},\boldsymbol{X}_{0:t-1}}\left(x_{t}|\left(\boldsymbol{n}_{t},j\right),\boldsymbol{x}_{0:t-1}\right)$ which requires $|\mathcal{S}^{\left(t-1\right)}|M$ multiplications.

Now consider the case where $n_{t,\ell}=1$ for some $\ell\in\mathcal{S}_{AR}$. Fix $t$ and $\ell$. Compute the sums $\sum\limits_{m=1}^{t}\alpha_{\boldsymbol{n}_{t}-\boldsymbol{1}+m\boldsymbol{e}_{\ell}}^{\left(t-1\right)}\left(\ell\right)$ for each $\boldsymbol{n}_{t}$ where $n_{t,\ell}=1$ and store them. After computing the sums, there are ${\left(|\mathcal{S}^{\left(t\right)}|-|\mathcal{S}^{\left(t-1\right)}|\right)}$ stored terms. Keep $t$ fixed, but allow $\ell$ to vary. Each of the stored sums is multiplied by $p_{\ell j}$ and $f^{\boldsymbol{\theta}}_{X_{t}|\boldsymbol{H}_{t},\boldsymbol{X}_{0:t-1}}\left(x_{t}|\left(\boldsymbol{n}_{t},j\right),\boldsymbol{x}_{0:t-1}\right)$ for $\ell\in\mathcal{S}_{AR}$ and $j\in\mathcal{S}$ which gives a total of $2Mk\left(|\mathcal{S}^{\left(t\right)}|-|\mathcal{S}^{\left(t-1\right)}|\right)$ multiplications. Thus, the total number of multiplications required is

[TABLE]

by the result $\binom{t}{m}\leq\frac{t^{m}}{m!}$. This can then be bounded by

[TABLE]

From which we see the complexity is bounded by $\mathcal{O}\left(M^{2}T^{k+1}{k^{k}}\right)$. ∎

To overcome possible underflow issues, we consider a normalised version of the algorithm. Define

[TABLE]

Lemma 7** (A normalised algorithm).**

Set $\widetilde{a}_{\boldsymbol{n}_{0}}^{\left(0\right)}={\alpha}_{\boldsymbol{n}_{0}}^{\left(0\right)}\left(j\right)$ from the simple algorithm. Then, for $t=1,\dots,T$ calculate

[TABLE]

Then the loglikelihood is given by

[TABLE]

Proof.

The definition of conditional densities gives

[TABLE]

where

[TABLE]

Using the same arguments as in the proof of Lemma 5, the right-hand side of (13) simplifies to

[TABLE]

$\text{when }n_{t,i}=1\text{ for some }i\in\mathcal{S}_{AR},$

[TABLE]

$\text{when }n_{t,\ell}\neq 1\text{ for all }\ell\in\mathcal{S}_{AR}.$ The summands in (14) and (15) can be written in the form

[TABLE]

which proves the result for the iterations. Equation (12) holds from the law of total probability. ∎

Lemma 8**.**

The complexity of the normalised forward algorithms, as given by the total number of multiplications, is less than $\mathcal{O}\left(M^{2}T^{k+1}k^{k}\right)$.

Proof.

The normalised algorithm is the same as the forward algorithm except that each term is divided by $\sum\limits_{\ell\in\mathcal{S}}\sum\limits_{\boldsymbol{n}\in\mathcal{S}^{\left(t-1\right)}}\widetilde{\alpha}_{\boldsymbol{n}}^{\left(t-1\right)}\left(\ell\right)$. The most efficient way to do this extra step is to do the division ${p_{ij}}/{\sum\limits_{\ell\in\mathcal{S}}\sum\limits_{\boldsymbol{n}\in\mathcal{S}^{\left(t-1\right)}}\widetilde{\alpha}_{\boldsymbol{n}}^{\left(t-1\right)}\left(\ell\right)}$ for $i,j\in\mathcal{S}$ first, which results in an additional $TM^{2}$ multiplications in total. ∎

3 A novel backward algorithm

The goal of the backward algorithm is to calculate the smoothed probabilities

[TABLE]

for $t=0,1,\dots,T$, $\boldsymbol{n}_{t}=\left(n_{1,t},\dots,n_{k,t}\right)\in\mathcal{S}^{\left(t\right)}$ and $i\in\mathcal{S}$. The smoothed probabilities are often of interest in their own right, but are also typically used to construct an EM algorithm. Recall that as a byproduct of the forward algorithm we obtain the filtered probabilities

[TABLE]

as well as the prediction probabilities

[TABLE]

for $j\in\mathcal{S}$, $\boldsymbol{n}_{t}\in\mathcal{S}^{\left(t\right)}$, and all $t=0,\dots,T$. These are the inputs to the backward algorithm.

Lemma 9** (A backward algorithm).**

The smoothed probabilities can be calculated using the following procedure. Set $\gamma_{\boldsymbol{n}_{T}}^{\left(T\right)}\left(i\right)=\widehat{\alpha}_{\boldsymbol{n}_{T}}^{\left(T\right)}\left(i\right),\text{ for all }i\in\mathcal{S},\text{ }\boldsymbol{n}_{T}\in\mathcal{S}^{\left(T\right)}.$ Then, for $t=T-1,T-2,\dots,0$ and for $\boldsymbol{n}_{t}\in\mathcal{S}^{(t)}$ calculate

[TABLE]

Proof.

Consider the event $\boldsymbol{H}_{t}=\left(\boldsymbol{n}_{t},i\right)$. By the definition of $\boldsymbol{n}_{t}$, when $i\in\mathcal{S}_{AR}$, then $\boldsymbol{n}_{t+1}=\boldsymbol{n}_{t}^{\left(-i\right)}+\boldsymbol{1}$, and when $i\in\mathcal{S}_{AR}^{c}$, then $\boldsymbol{n}_{t+1}=\boldsymbol{n}_{t}+\boldsymbol{1}$. Thus, when $\boldsymbol{H}_{t}=\left(\boldsymbol{n}_{t},i\right)$ is known, then $\boldsymbol{n}_{t+1}$ is also known. As a result,

[TABLE]

for $t=0,1,\dots,T-1$, $\boldsymbol{n}_{t}\in\mathcal{S}^{\left(t\right)}$ and $i\in\mathcal{S}$. Since the following arguments are the same for both cases, $i\in\mathcal{S}_{AR}$ and $i\in\mathcal{S}_{AR}^{c}$, for notational convenience, let $\boldsymbol{n}$ take the value $\boldsymbol{n}_{t}^{\left(-i\right)}+\boldsymbol{1}$ when $i\in\mathcal{S}_{AR}$ and the value $\boldsymbol{n}_{t}+\boldsymbol{1}$ when $i\in\mathcal{S}_{AR}^{c}$. The summands on the right hand side of (16) can be written as

[TABLE]

where the last equality holds since $\boldsymbol{H}_{t+1}$ is independent of $\boldsymbol{x}_{0:t}$ given $\boldsymbol{H}_{t}$. Now, noting that $\boldsymbol{x}_{t+1:T}$ is independent of $\boldsymbol{H}_{t}$ given $\boldsymbol{H}_{t+1}$ and $\boldsymbol{x}_{0:t}$, then the right-hand side of (17) equals

[TABLE]

Writing out $\boldsymbol{n}$ explicitly for the two cases completes the proof. ∎

Lemma 10**.**

The total complexity of the backward algorithm in Lemma 9, as measured by the total number of multiplications, is less than $\mathcal{O}\left(M^{2}T^{k+1}k^{k}\right).$

Proof.

First, for each $t\in\{T-1,...,0\}$ we need to calculate the ratio

[TABLE]

for every corresponding $\boldsymbol{n}_{t}\in\mathcal{S}^{\left(t\right)}$ and $j\in\mathcal{S}$. This costs $M|\mathcal{S}^{\left(t\right)}|$ multiplications. This quantity is independent of $i$, thus only needs to be done once for a given $t$ if we save the results.

Now consider $t$, $i$ and $\boldsymbol{n}_{t}$ fixed. The multiplication of $p_{ij}$ and

[TABLE]

is done for every $j\in\mathcal{S}$ which costs $M$ multiplications. The sum over $j\in\mathcal{S}$ results in a single term, which is then multiplied by the corresponding $\widehat{\alpha}_{\boldsymbol{n}_{t}}^{\left(t\right)}\left(i\right)$, and this costs an additional 1 multiplication. We do this for all $i\in\mathcal{S}$ and $\boldsymbol{n}_{t}\in\mathcal{S}^{\left(t\right)}$, which costs $\left(M+1\right)M|\mathcal{S}^{\left(t\right)}|$ multiplications.

So, for a given $t$ we execute $M|\mathcal{S}^{\left(t\right)}|+\left(M+1\right)M|\mathcal{S}^{\left(t\right)}|$ multiplications. This is done for every $t=0,\dots,T-1$, so the total number of multiplications is

[TABLE]

where we have used similar arguments to Lemma 6 to bound the complexity. ∎

Of importance to the next section, note that we can obtain from the backward algorithm, the smoothed probabilities

[TABLE]

4 A novel EM algorithm

Here we show how the output from the backward algorithm can be used to implement an exact, computationally feasible EM algorithm for MRS models with independent regimes.

4.1 The E-step

Recall that the EM algorithm [7] is an iterative procedure, alternating between an expectation step and a maximisation step. In the expectation step the function $Q\left(\boldsymbol{\theta},\boldsymbol{\theta}_{n}\right)$ is constructed as

[TABLE]

where $\boldsymbol{R}=\left(R_{0},\dots,R_{T}\right)$ is a sequence of the hidden Markov chain $\{R_{t}\}$, and $\left(\boldsymbol{H}_{0},\dots,\boldsymbol{H}_{T}\right)$ is a sequence of the corresponding augmented hidden process $\{\boldsymbol{H}_{t}\}$. The information contained in the sequences $\boldsymbol{R}$ and $\left(\boldsymbol{H}_{0},\dots,\boldsymbol{H}_{T}\right)$ is entirely equivalent, but we opt for the latter representation to remain consistent with, and emphasise the place of, the work in the previous sections. In the M-step of the algorithm, the maximisers $\operatorname*{arg\,max}\limits_{\boldsymbol{\theta}\in\Theta}Q\left(\boldsymbol{\theta},\boldsymbol{\theta}_{n}\right)$ are found.

For MRS models $Q\left(\boldsymbol{\theta},\boldsymbol{\theta}_{n}\right)$ can be written as

[TABLE]

Using the augmented hidden Markov chain, $\{\boldsymbol{H}_{t}\}_{t\in\mathbb{N}}$, (19) can be written in such a way that the function $Q$ is computationally feasible. First note that, given $\boldsymbol{H}_{t}$ and $\boldsymbol{X}_{0:t-1}$, $X_{t}$ is independent of $\boldsymbol{H}_{s}$ for $s\neq t$ which allows the function $\log f^{\boldsymbol{\theta}}_{\boldsymbol{X}_{0:T}|\boldsymbol{H}_{0},\dots,\boldsymbol{H}_{T}}\left(\boldsymbol{x}_{0:T}|\boldsymbol{H}_{0},\dots,\boldsymbol{H}_{T}\right)$ to be written as

[TABLE]

Since $f^{\boldsymbol{\theta}}_{X_{t}|\boldsymbol{H}_{t},\boldsymbol{X}_{0:t-1}}\left(x_{t}|\left(\boldsymbol{n}_{t},j\right),\boldsymbol{x}_{0:t-1}\right)=f^{\boldsymbol{\theta}}_{X_{t}|N_{t,j},R_{t},\boldsymbol{X}_{0:t-1}}\left(x_{t}|n_{t,j},j,\boldsymbol{x}_{0:t-1}\right)$, and similarly for $f^{\boldsymbol{\theta}}_{X_{0}|\boldsymbol{H_{0}}}\left(x_{0}|\left(\boldsymbol{n}_{t,j},j\right)\right)$, the expression (21) simplifies to

[TABLE]

Taking the expectation of (23) with respect to the distribution $f^{\boldsymbol{\theta}_{n}}_{\boldsymbol{H}_{0},\dots,\boldsymbol{H}_{T}|\boldsymbol{X}_{0:T}}$ (equivalently the distribution $f^{\boldsymbol{\theta}_{n}}_{\boldsymbol{R}|\boldsymbol{X}_{0:T}}$) gives

[TABLE]

Using similar arguments, $\mathbb{E}\left[\log\mathbb{P}^{\boldsymbol{\theta}}\left(\boldsymbol{H}_{0},\dots,\boldsymbol{H}_{T}\right)|\boldsymbol{x}_{0:T};\boldsymbol{\theta}_{n}\right]$ is found to be

[TABLE]

where $\eta_{ij}$ is the random variable counting the number of transitions from state $R_{t-1}=i$ to state $R_{t}=j$ in the sequence $\boldsymbol{R}$. The expectation $\mathbb{E}\left[\eta_{ij}|\boldsymbol{x}_{0:T};\boldsymbol{\theta}_{n}\right]$ can be calculated as

[TABLE]

So, the function $Q$ is

[TABLE]

Lemma 11**.**

The joint probabilities are given by

[TABLE]

This proof follows similar arguments to those in [20], which develops algorithms for MRS models with dependent regimes.

Proof.

For the case $i\in\mathcal{S}_{AR}$, note that $N_{t,i}=1$ if and only if $R_{t-1}=i$ and we are done.

When $i\in\mathcal{S}_{AR}^{c}$ all counters in $\boldsymbol{n}_{t}$ are different from 1, so $\boldsymbol{n}_{t}-\boldsymbol{1}\in\mathcal{S}^{\left(t-1\right)}$. Thus

[TABLE]

The last equality holds since, given $R_{t}$ and $\boldsymbol{N}_{t}$, then $\boldsymbol{x}_{t:T}$ is independent of $R_{t-1}$. Focusing on the right-most term in Equation (26),

[TABLE]

where $p_{ij}^{\left(n\right)}$ is the parameter $p_{ij}$ in $\boldsymbol{\theta}_{n}$; the second equality holds since, given $R_{t-1}$, then $R_{t}$ is independent of $\boldsymbol{N}_{t}$ and $\boldsymbol{X}_{0:t-1}$.

Now, notice that

[TABLE]

since $i\in\mathcal{S}_{AR}^{c}$, and that

[TABLE]

with the sum in the denominator being over $k\in\mathcal{S}_{AR}^{c}$ since only when $k\in\mathcal{S}_{AR}^{c}$ is $\boldsymbol{n}_{t}-\boldsymbol{1}\in\mathcal{S}^{(t-1)}$ defined; this completes the proof. ∎

4.2 The M-step

Next, the maximisers, $\boldsymbol{\theta}_{n+1}=\operatorname*{arg\,max}\limits_{\boldsymbol{\theta}\in\Theta}Q\left(\boldsymbol{\theta},\boldsymbol{\theta}_{n}\right)$, are needed. The maximisers for the parameters of each regime are generally problem specific, but the maximisers for the parameters $p_{ij}$, $i,j\in\mathcal{S}$, can be derived in general. By the work of [11],

[TABLE]

However, note that to get this analytic update for the $p_{ij}^{\left(n+1\right)}$ parameters, terms involving $\pi_{j}$ in Equation (24) have been treated as if they are unrelated to $p_{ij}$, $i,j\in\mathcal{S}$. However, this is not true when $\boldsymbol{\pi}$ is specified as the stationary distribution of the process $\{R_{t}\}$, but holds for other cases, such as when $\boldsymbol{\pi}$ is some predetermined distribution, or when $\boldsymbol{\pi}$ is specified as a parameter to be inferred. Nonetheless, this simplification is appropriate if we assume that, as the sample size grows, the contribution of terms involving $R_{0}$ become insignificant.

4.3 Model-specific M-step updates

In electricity price models, it is common to specify spike or drop regimes as either shifted-Gamma, shifted-log-normal, or occasionally a Gaussian distribution. Here we derive M-step updates for these regimes.

Corollary 12**.**

Suppose Regime $i$ is i.i.d. $N\left(\mu_{i},\sigma_{i}^{2}\right)$. The M-step updates $\mu_{i}^{\left(n+1\right)}$ and $\left(\sigma^{\left(n+1\right)}_{i}\right)^{2}$, for $n\geq 0$, are

[TABLE]

Proof.

The results holds after differentiating the $Q\left(\cdot,\boldsymbol{\theta}_{n}\right)$ and solving for zeros. That $\mu_{i}^{\left(n+1\right)}$ is a maximiser is shown by the second derivative test. Furthermore, $\left(\sigma_{i}^{\left(n+1\right)}\right)^{2}$ can be shown to be a maximiser by comparing the value of $Q\left(\boldsymbol{\theta},\boldsymbol{\theta}_{n}\right)$ when $\sigma_{i}^{2}=\left(\sigma_{i}^{\left(n+1\right)}\right)^{2}$ to the value of $Q\left(\boldsymbol{\theta},\boldsymbol{\theta}_{n}\right)$ evaluated at any other value of $\sigma_{i}^{2}$, and utilising the inequality $1-{1}/{y}\leq\log y$. ∎

At time $t$, if $R_{t}=i$ is a shifted-Gamma regime so that $X_{t}-q_{i}\sim Gam\left(\mu_{i},\sigma_{i}^{2}\right)$, where $q_{i}$ is known, the M-step is not completely analytic. However, the dimension of the maximisation problem can be reduced from 2-dimensional to 1-dimensional via the following corollary.

Corollary 13**.**

Suppose Regime $i$ follows an i.i.d. shifted-Gamma distribution, that is, if $X_{t}$ is from Regime $i$, then $\left(X_{t}-q_{i}\right)\sim Gam\left(\mu_{i},\sigma_{i}^{2}\right)$, and suppose the parameter $q_{i}$ is known. The M-step update for the scale parameter, $\left(\sigma_{i}^{\left(n+1\right)}\right)^{2}$, as a function of $\mu_{i}^{\left(n+1\right)}$, is

[TABLE]

The update for $\mu_{i}$ is then found by finding

[TABLE]

where $\Gamma\left(\cdot\right)$ is the Gamma function.

Proof.

The result follows after differentiating $Q$ with respect to $\sigma_{i}^{2}$, and solving for the stationary point, which is a maximum by the second derivative test. ∎

Corollary 14**.**

Suppose Regime $i$ follows i.i.d. shifted-log-normal dynamics, that is, if $X_{t}$ is from Regime $i$, then $\log\left(X_{t}-q_{i}\right)\sim N\left(\mu_{i},\sigma_{i}^{2}\right)$, and suppose the parameter $q_{i}$ is known. The M-step updates are

[TABLE]

Proof.

The proof is similar to the proof of Corollary 12. ∎

Note that Corollaries 13 and 14 assume the parameter $q_{i}$ is known. This is necessary for the shifted-log-normal distribution [12] and the shifted-Gamma distribution when the shape parameter $\mu_{i}$ is less than 1 [19]. Furthermore, for the shifted-Gamma distribution, [19] observe that related issues arise when $\mu$ is near 1, and advise against maximum likelihood estimation of $q_{i}$ when $\mu<2.5$. Simulations suggest this is also good advice when fitting MRS models with shifted-Gamma regimes [22].

In electricity price modelling literature it is common to specify a ‘base regime’ as an AR(1) process. In existing literature the AR(1) regimes are assumed to evolve at every time $t$ but are only observed when in that regime. Another possibility is that the AR(1) processes evolve only when they are observed, that is, define $\tau_{i}\left(t\right)=\sum\limits_{\ell=0}^{t}\mathbb{I}\left(R_{\ell}=i\right)$, and the AR(1) process in Regime $i$ as $\{B_{\tau_{i}\left(t\right)}\}_{\tau_{i}\left(t\right)\in\mathbb{N}}$. Our algorithms are applicable to both specifications; however, for simplicity, here we treat the former only. For more details on the latter specification, see [22].

Corollary 15**.**

If Regime $i$ is an AR(1) regime of an MRS model, the M-step of the EM algorithm can be executed as follows. The updates $\alpha_{i}^{\left(n+1\right)}$ and $\left(\sigma_{i}^{\left(n+1\right)}\right)^{2}$ as functions of $\phi_{i}^{\left(n+1\right)}$ are

[TABLE]

The M-step update for $\phi_{i}^{\left(n+1\right)}$ is given by

[TABLE]

where

[TABLE]

Proof.

Differentiate $Q$ with respect to $\alpha_{i}$ and $\sigma_{i}^{2}$ and solve for when the derivative is zero. The parameter $\alpha_{i}$ can be shown to be a maximiser by the second derivative test, and $\sigma_{i}^{2}$ can be seen to be a maximiser using the same argument as used in the proof of Corollary 12. Next, substitute the maximisers, $\alpha_{i}^{\left(n+1\right)}\left(\phi_{i}^{\left(n+1\right)}\right)$ and $\sigma_{i}^{\left(n+1\right)}\left(\phi_{i}^{\left(n+1\right)}\right)^{2}$, into (24) and collect all terms involving $\phi_{i}$, to give the function $g$. That we need to search for the global maximiser of $g$ on the interval $\left(-1,1\right)$ only comes from the fact that we have assumed Regime $i$ is a stationary or mean-reverting process, in which case $|\phi_{i}|<1$ is a necessary condition. ∎

Remark 4.1: Truncation

The forward and backward algorithms can be computationally costly when $T$ and/or $k$ are large, so it may be preferential (or necessary) to truncate the problem. We suggest that the memory of each of the AR(1) processes be truncated. That is, for all $\boldsymbol{n}_{t}$ such that $n_{t,i}>D-1$ for some $i\in\mathcal{S}_{AR}$, we let $f_{X_{t}|\boldsymbol{H}_{t},\boldsymbol{X}_{0:t-1}}\left(x_{t}|\left(\boldsymbol{n}_{t},i\right),\boldsymbol{x}_{0:t-1}\right)=\widetilde{f}_{X_{t}|R_{t},\boldsymbol{X}_{0:t-1}}\left(x_{t}|i,\boldsymbol{x}_{0:t-1}\right),$ where $\widetilde{f}_{X_{t}|R_{t},\boldsymbol{X}_{0:t-1}}\left(x_{t}|i,\boldsymbol{x}_{0:t-1}\right)$ is a density that does not depend on $\boldsymbol{n}_{t}$. An appropriate choice of $\widetilde{f}$ will be problem-specific. This truncation is equivalent to truncating the state space of $\boldsymbol{H}_{t}$ so that $N_{t,i}\in\{1,\dots,D\}$, and adjusting the transitions of $\boldsymbol{H}_{t}$ so that the counters remain at $D$, rather than continuing to increase as they would in the original process. It can be shown that the complexities of the truncated algorithms are $\mathcal{O}(M^{2}D^{k}Tk^{k})$. For processes that decay to stationary between times at which they are observed (such as the AR(1) processes in independent-regime MRS models used in this paper, see Section 4.3), then taking $\widetilde{f}_{X_{t}|R_{t},\boldsymbol{X}_{0:t-1}}\left(x_{t}|i,\boldsymbol{x}_{0:t-1}\right)$ as the stationary distribution in Regime $i$ is a logical choice. Furthermore, $D$ should be chosen large enough so that there is a low probability that $\{R_{t}\}$ ever visits any specific state for more than $D$ consecutive transitions and of course $D>k$.

5 Simulation studies

We perform a simulation study to examine the properties of the algorithms. The models used in the simulations are the following:

[TABLE]

where $\{B_{t}\}$ is an AR(1) process defined by $B_{t}=0+0.75B_{{}_{t}-1}+\varepsilon_{t}$, $\varepsilon_{t}\sim$ i.i.d. $N(0,1)$, $S_{t}\sim$ i.i.d N$\left(0,1\right)$ and $\{R_{t}\}$ is a Markov chain with state space $\mathcal{S}=\{1,2\}$, transition matrix entries $p_{11}=p_{22}=0.9$, and initial distribution $\left(0.5,0.5\right)$;

[TABLE]

where $\{B_{t}^{\left(i\right)}\}$ are independent AR(1) process defined by $B_{t}^{\left(1\right)}=0+0.9B_{t-1}^{\left(1\right)}+\varepsilon^{\left(1\right)}_{t}$ and $B_{t}^{\left(2\right)}=0+0.4B_{t-1}^{\left(2\right)}+\varepsilon_{t}^{\left(2\right)}$, $\varepsilon_{t}^{\left(i\right)}\sim$ i.i.d. $N(0,1)$ $i=1,2$, and $\{R_{t}\}$ is a Markov chain with state space $\mathcal{S}=\{1,2\}$, transition matrix entries $p_{11}=p_{22}=0.6$, and initial distribution $\left(0.5,0.5\right)$.

To investigate the bias and consistency of the MLE, 20 independent realisations of Models 1 and 2 were simulated for $T=50,100,200$ and $400$, and the EM algorithm used to find the MLE. The terminating criteria for the algorithms was to stop when either the increase in the likelihood, or the step-size, as measure by $|\boldsymbol{\theta}_{n+1}-\boldsymbol{\theta}_{n}|_{\infty}$ was less that $1.5\times 10^{-}8$. To attempt to avoid local maxima, the EM algorithm was initialised from 100 randomised values centred around the true parameters. Of the corresponding 100 terminating points of the EM algorithm, the parameters that achieved the highest loglikelihood value were kept. Figure 3 shows box plots of the 20 terminating points of the EM algorithm, one for each simulation, as well as the log-runtime and value of the loglikelihood. For both models the MLE appears to be converging to the true parameter value as sample sizes increase. Generally, there appears to be a much larger variation in the MLEs for Model 2 than for Model 1. This could be because the inference problem for Model 2 is harder, as the regimes in Model 2 are more similar than they are in Model 1, or because of the nature of the hidden Markov chain is such that there is a lower probability of remaining in each regime ($p_{11}=p_{22}=0.6$ for Model 2, compared to $p_{11}=p_{22}=0.9$ for Model 1), or both. Regarding the run time, for data sets of length 100 and greater, the empirical results suggest the complexity is approximately $\mathcal{O}\left(T^{1.02}\right)$ for Model 1 and $\mathcal{O}\left(T^{1.70}\right)$ for Model 2, which agrees with, and is significantly better than, our theoretical upper bound.

To investigate the truncation method, we used the same data sets simulated above for $T=400$ and applied our truncated algorithm to this data with various truncation levels, $D=5,10,20,40$. The termination criteria for the EM algorithm was to stop when either the increase in the likelihood, or the step-size, as measure by $|\boldsymbol{\theta}_{n+1}-\boldsymbol{\theta}_{n}|_{\infty}$ was less that $1.5\times 10^{-8}$. Figure 4 plots the $\log_{10}$ of the absolute value of the difference between the parameters recovered by the truncated algorithm and the full algorithm. Figure 4 also shows $\log_{10}$ of the absolute value of the difference in the loglikelihoods achieved by the truncated algorithm and the full algorithm. Figure 4 also plots of the relative runtime of the truncated algorithm compared to the full algorithm, that is, the ratio of the time taken for the truncated algorithm compared to the full algorithm. The run times for the truncated algorithms are significantly lower than the full algorithm. For these data sets the truncated algorithm performs reasonably well, even when limiting the memory of the counters to just 5 time steps. For a truncation level of 40, the errors are of the order $10^{-6}$ to $10^{-10}$ for both models, which is similar to the stopping criteria for the EM algorithm which is of the order $10^{-8}$. For Model 1 the truncation appears to have a more significant negative effect on the parameter estimates compared to Model 2. This could be because the probability of staying in each regime is higher in Model 1 than in Model 2, and therefore the probability of remaining in one regime for more that $D$ consecutive transitions is lower in Model 2. It is likely that the error due to truncation of the algorithm is insignificant compared to statistical error of parameter estimates.

To investigate the convergence properties of the EM algorithm we used five of the simulated data from above with $T=400$ and used the truncated EM algorithm, with a truncation level of $D=40$ to search for maxima. For each of the five simulated datasets we ran the EM algorithm 50 times, sampling initial values for the algorithm independently each time and used the same terminating criteria as before. The sampling distributions for the initial parameter values of the algorithm are summarised in Table 1. For Model 2 we have ordered the terminating values of the EM algorithm such that $\phi_{1}>\phi_{2}$ for identifiability.

Observing Figure 5 we see that, for most of the simulations and most of the starting values, the EM algorithm finds a single maxima. In Figure 5 the black cross represents the point which achieved the highest maximum. We refer to this point as the optimal parameter value. The proportion of times that the algorithm converged to the optimal parameter value is reported in Table 2. Simulation 4 of Model 1 is outstanding since only 14% of initial values resulted in the algorithm finding the optimal value. In Figure 5 there are some instances when the algorithm terminates at a suboptimal point; in particular simulated datasets 3 and 4 for Model 1. As shown in. Figure 5, for Model 1 the optimal parameter set appears to reasonably estimate the true parameters for all simulated dataset.

Estimating Model 2 is much harder since the regimes are both very similar. Figure 5 shows that the optimal parameter set estimates the true parameters reasonably for simulated datasets 2-5, but not for simulated dataset 1. The behaviour of the algorithm for simulated dataset 1 for Model 2 is particularly interesting. For this dataset, the algorithm terminates at one of two distinct locations. One of these terminating points is at $p_{11}=1$ and hence Regime 1 is absorbing. This means the algorithm has converged to a point where all but the first observation are captured by Regime 1. As $p_{11}$ tends to $1$ the algorithm is able to send $\sigma_{2}^{2}\to 0$ to achieve arbitrarily large values of the likelihood and the model is unidentifiable. To prevent this behaviour, we suggest that the parameters $\sigma_{i}^{2}$ and $p_{ij}$ are restricted so that they are away from the boundary. Indeed, for simulated dataset 1 of Model 2, if we take the terminating values of the algorithm which lie away from the boundary, then we get reasonable estimates of the true parameters. For more discussion see [22].

6 An application to South Australian wholesale electricity market

The dataset consists of 81,792 half-hourly spot prices from the South Australian electricity market (available at the AEMO website [3]) for the period 00:00 hours, $1^{\text{st}}$ of January 2013, to 23:30 hours, $31^{\text{st}}$ of September 2017. Note that this dataset contains a period of 14 days over which the market was suspended from 4:00pm, on the $28^{\text{th}}$ of September until 10:30pm on the $11^{\text{th}}$ of October. During this period prices were set by AEMO. We ignore this fact in our modelling and include them in the data set anyway.

Following a common practice in the literature we model daily average prices and thus we have a dataset of 1,704 daily average price observations to which we fit our model. The data that we model is plotted in Figure 6. To model the South Australian wholesale electricity market, we break the price process up in to two components, $P_{t}=S_{t}+X_{t}$, where $P_{t}$ is the price on day $t$, $S_{t}$ is a deterministic trend component, and $X_{t}$ is a stochastic component which is to be modelled by an independent regime MRS model.

Electricity spot prices exhibit seasonality on daily, weekly, and longer scales. To capture this multi-scale seasonality, the trend component consists of two parts: a short-term component, $g_{t}$, and a long-term component, $h_{t}$, so $S_{t}=g_{t}+h_{t}$. We model the long-term component, $h_{t}$, using wavelet filtering since it has been shown to perform well for this application [18], and use Daubechies 24 wavelets and a level 6 approximation [18]. We use the short-term component, $g_{t}$, to capture the mean price for different days of the week and indicator functions to model this:

[TABLE]

where $\beta_{\text{Mon}}$, $\beta_{\text{Tue}},$ …, $\beta_{\text{Sun}}$ are the mean deviations from the long-term trend price on Monday, Tuesday, …, Sunday, respectively.

Following [18] we use the RFP (recursive filter on prices) method to estimate the seasonal component in the presence of extreme observations. The method first uses the raw price series to estimate the trend model, then removes this from the data. Next, the standard deviation of these altered prices is estimated, and any prices that are more than three standard deviations from the current estimate of the trend are replaced with the value of the estimated trend at that point. The procedure then re-estimates the trend component on the original data set with spikes removed.

We consider the following four models for the stochastic component:

[TABLE]

where $B_{t}=\alpha+\phi B_{t-1}+\sigma_{1}\varepsilon_{t}$ is an AR(1) process with $\varepsilon_{t}\sim$ i.i.d. N$\left(0,1\right)$, and $S_{t}-q_{3}$, where $q_{3}$ is a shifting parameter, follows either a Gamma distribution (M1-Gamma), or a log-normal distribution with parameters $\mu_{2}$ and $\sigma_{2}^{2}$ (M1-LN). The last two models introduce a ‘drop’ regime as well:

[TABLE]

where $B_{t}$ and $S_{t}$ are as above, and $-D_{t}+q_{1}$, where $q_{1}$ is a shifting parameter, follows a log-normal distribution with parameters $\mu_{3}$ and $\sigma_{3}^{2}$. The use of the shifting parameters $q_{1}$ and $q_{3}$ was proposed by [15]. As previosuly mentioned, estimating the shifting parameters of these distribution is known to be a difficult task [19], so we fix $q_{1}$ and $q_{3}$ as the first and third quantiles of the detrended data, as suggested by [17].

The models were fitted to the detrended data using the truncated EM algorithm with truncation level $D=56$ (8 weeks). The Bayesian Information Criterion (BIC) values of these models are reported in Table 3. From the BIC values we choose model M1-LN, with an AR(1) base regime and a log-normal spike regime. The parameters for this model are reported in Table 4. Of course, for a rigorous treatment of this modelling problem, model assumptions should be checked and we should not rely solely on the BIC. Using the smoothed probabilities obtained while fitting model M1-LN, the prices can be classified into which regime is most likely. This is shown in Figure 6, where prices are highlighted in red if $\mathbb{P}^{\widehat{\boldsymbol{\theta}}}\left(R_{t}=2|\boldsymbol{x}_{0:T}\right)>0.5$, where $\widehat{\boldsymbol{\theta}}$ is the MLE.

7 Conclusions

In this paper we have developed novel techniques for independent-regime MRS models. Specifically, we consider models that are a collection of independent AR(1) processes, where only one process is observed at each time $t$, and which regime is observed is determined by a hidden Markov chain. We develop forward, backward and EM algorithms for these models, and show that the methods we develop here can outperform the existing method of inference used in the electricity price modelling literature, the EM-like algorithm [17].

The construction of these methods relies on the idea of augmenting the hidden Markov chain with a set of counters, which keep track of the number of transitions since the last visit to AR(1) regimes. The forward algorithm can be used to evaluate the likelihood and filtered and prediction probabilities, which are used as inputs to the backward algorithm. The backward algorithm is used to evaluate the smoothed probabilities. Together, the forward-backward procedure executes the E-step of the EM algorithm. We showed that the complexity of the forward and backward algorithms is $\mathcal{O}\left(M^{2}T^{k+1}k^{k}\right)$, where $M$ is the number of regimes in the model, $T$ the length of the observed sequence, and $k$ the number of AR(1) regimes in the model. This complexity may be impractically large if $T$ or $k$ are large, so we introduce an approximation where the memory of AR(1) processes in the model is truncated. These truncated methods are $\mathcal{O}\left(M^{2}D^{k}Tk^{k}\right)$, where $D$ is the memory of the counters.

Simulation suggests that the MLE found via the EM algorithm is consistent, and that, even for models with regimes with similar characteristics (such as Model 2), the MLE is still a reasonable estimator when the sample size is large enough (400 observations in this case). Simulations also suggest that the truncation method is a reasonable approximation to the full likelihood, while improving runtime and memory requirements significantly. As is typical for hill-climbing algorithms, there is the possibility of the algorithm terminating at sub-optimal values depending on the initial values of the algorithm. We explored some of this behaviour via a simulation study and found, at most, two local maxima of the likelihood function. Of course, the behaviour of the algorithm is going to be model- and data-specific. This simulation study also highlighted possible identifiability issues which can arise. However, these can be rectified by restricting parameters away from the boundary.

Lastly, we apply our methods to estimate four models for the South Australian wholesale electricity market and find that a 2-regime model with an AR(1) base regime and shifted log-normal spike regime is best as measured by the BIC. We also demonstrate how prices can be classified into regimes using the smoothed probabilities.