6-10-2025 - MS-AR Implementation

Markov Switching Auto Regression (MS-AR) Algorithm

Since the last update I’ve spent a few late nights finishing off the full code implementation of the Markov Switching Auto Regression algorithm (MS-AR), this means not just finding the $\mu$ and $\sigma$ values for each of the states but fitting a regression model with a pre-specified auto regressive order to each state, so each state develops a unique relationship between periodic returns within the data.

This implementation is built on top of the Expectation Maximisation (EM) algorithm which is an algorithm that iteratively tunes each states parameters to converge to $n$ $\theta$ weights for each states’ regression function. So as a result the algorithm is essentially built in 3 different sections:

---

Initialisation

Due to the nature of the EM algorithm its circular so to begin we estimate some starting values that can then be used to kickstart the iterative process of tuning parameters - think of this as the starter motor in a car, it gets the engine going initially. So for this we just need to estimate initial values for each states regressive function, this can be done through a standard clustering algorithm where we cluster the data into $n$ clusters. Next we perform Ordinary Least Squares regression (OLS) on each cluster individually, this is like performing Weighted Least Square regression but because we have used hard probabilities - ie a value either is or isn’t in the given cluster - theres no need to multiply by probabilities of being in each state/cluster so it’s simpler just applying the OLS algorithm to each cluster individually due to having the same effect.

The formula for OLS regression is:

${\hat{\mathbf{\beta }}}_{\text{OLS}}=({\mathbf{X}}^{\text{T}}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\text{T}}\mathbf{y}$

Now we have estimated parameters for each state’s regression algorithm we can produce a standard deviation, for this we calculate a prediction for each periodic return based on the each state’s regression algorithm - so producing a matrix with a shape of $n\times m$ where $n$ is the number of periodic returns in our data set and $m$ is the number of identified states - then for every value in each state we calculate the residual of the prediction then finally we filter each residual based on what cluster its target value is and use all residuals belonging to target values in cluster $n$ to calculate the initial standard deviation for state $n$ hence giving us the initial values we need for for the EM algorithm - regression function coefficients ($\theta$ values) and standard deviation ($\sigma$) for each state

---

Expectation Step

The aim of the expectation step is to get the probability that each value is in each state - so the output would be an $n\times m$ matrix where $n$ is the number of individual periodic return data points - the size of the autoregressive order (ARO) and $m$ is the number of states being used within the system.

First we generate a new regressive prediction for each state using each state’s regressive coefficients for every data value in the time series data. This is achieved by using a feature matrix - $\mathbf{X}$ - that has dimension $([\text{number of periodic return data points - ARO}]\times \text{ARO}+1)$. Each row within this matrix has the structure $\{1.0,y(t-1),y(t-2),...,y(t-\text{ARO})\}$, where $y$ is the target data value - in this case periodic returns - which leaves us with a matrix of time lagged data. By using this matrix it makes it much simpler to calculate the prediction for each state based on the state’s regressive function coefficients as we just perform a dot product of the time lagged data - $\mathbf{X}$ - by the vector of the state’s regression function coefficients.

After we have these predictions for each value in every state we then calculate the residual of the prediction in each state to the actual target value - $y$. Then using the residuals we can use the current state’s calculated $\sigma$ value along with an assumed $\mu$ of 0 and plug them into a continuous data probability density function to calculate the likelihood that the given residual is in the given state. The formula for the probability density function (PDF) for a normal distribution is:

$p(x)=\frac{1}{\sigma \sqrt{2\pi }}{e}^{-\frac{1}{2}{\left(\frac{x-\mu }{\sigma }\right)}^2}$

The final step is to then row normalise the likelihoods across each state and scale them to sum to 1 whilst maintaining their relative difference hence giving us the probability each residual is in each state and therefore rounding out the expectation step of the algorithm.

---

Maximisation Step

This is the step where we re-tune the regression function coefficients and weighted standard deviation based on the probabilities produced in the expectation step in the algorithm. First we take the the probability values generated by the expectation step for each residual in every state. Then we calculate the new $\theta$ values for the given state by using a Weighted Least Squares (WLS) function that takes the Weights (probability each residual is in the given state), $\mathbf{W}$, the time lagged periodic returns matrix from the expectation step - $\mathbf{X}$, and finally the target/actual value for each data point in the time series data - $\mathbf{y}$. Then instead of having to iteratively tune the regression function coefficients we just use a series of matrix multiplications to instantly calculate the optimal regression coefficients.

The function for WLS:

$\hat{\mathbf{\beta }}=({\mathbf{X}}^{\text{T}}\mathbf{WX}{)}^{-1}{\mathbf{X}}^{\text{T}}\mathbf{Wy}$

Then $\hat{\mathbf{\beta }}$ represents the $\theta$ values for the regression function of the given state which we use to calculate new residuals for every value in each state, then the final step is to take these new residuals, multiply them by their corresponding weights calculated in the expectation step and then find the standard deviation of the resulting values for each state to get the given state’s weighted standard deviation. The maximisation step then returns the new $\theta$ values and standard deviation for each state.

The code for this can be found on my git hub page within the market-trading-algorithm repository.