A Hidden Markov Model (HMM) is a statistical model used to represent systems where observed data are generated by an underlying, unobserved (hidden) sequence of states following a Markov process. It is widely applied in fields such as speech recognition, biology, finance, and economics.

The key concepts related to HMMs are listed in the table below. We will explain each concept as we walk through a simplified stock market example.

Our simplified stock market example begins by specifying the hidden states (S), observable outputs, and the observation sequence (O), as follows:

• 2 hidden states:
  • Bull market
  • Bear market
• 2 observable outputs:
  • Up day (stock price rises)
  • Down day (stock price falls)
• Observation sequence (based on 3 days of market movement):
  • Day 1: up
  • Day 2: down
  • Day 3: up
    (So: o=[up, down, up])

The model parameters are initialized as follows (with explanations provided below the data matrices):

The initial probabilities above assume a 60% chance that the market is in a Bull state and a 40% chance that it is in a Bear state.

The transition matrix indicates that the probability of a Bull market remaining Bull is 80%, transitioning from Bull to Bear is 20%, from Bear to Bull is 30%, and from Bear to Bear is 70%.

The emission matrix shows that, given the market is in a Bull state, the probability of observing “Up” is 90% and “Down” is 10%. If the market is in a Bear state, the probability of observing “Up” is 20% and “Down” is 80%.

The implementation of an HMM involves three key types of probabilities: forward, backward, and posterior (or state) probabilities. Forward probabilities represent the probability of observing the sequence up to time t and being in state i at time t. Backward probabilities represent the probability of observing the future observations from time t+1 onward, given state i at time t. Posterior probabilities combine both to estimate the probability of being in state i at time t, given the entire sequence of observations.

First, we calculate the forward probabilities (α).

Next, we calculate the backward probabilities (β).

Finally, we calculate the posterior probabilities.

The transition matrix, emission matrix, and initial state probabilities are iteratively updated during training using the Baum-Welch algorithm. Once the model is trained, these parameters are fixed and used for inference tasks such as decoding or state classification.

For brevity, we do not show the detailed update calculations here, but the accompanying Python code includes the full parameter update logic.

The example above can be implemented in Python as follows:

import numpy as np
import pandas as pd

# Observations: Up=0, Down=1, Up=0
observations = np.array([0, 1, 0])
T = len(observations)
N = 2  # Hidden states: Bull, Bear
M = 2  # Observation symbols: Up, Down

# Initial model parameters
pi = np.array([0.6, 0.4])
A = np.array([
    [0.8, 0.2],  # Bull to Bull, Bull to Bear
    [0.3, 0.7]   # Bear to Bull, Bear to Bear
])
B = np.array([
    [0.9, 0.1],  # Bull emits Up, Down
    [0.2, 0.8]   # Bear emits Up, Down
])

# Step 1: Forward algorithm (alpha)
alpha = np.zeros((T, N))
alpha[0] = pi * B[:, observations[0]]
for t in range(1, T):
    for j in range(N):
        alpha[t, j] = np.sum(alpha[t-1] * A[:, j]) * B[j, observations[t]]

# Step 2: Backward algorithm (beta)
beta = np.zeros((T, N))
beta[T-1] = 1
for t in range(T-2, -1, -1):
    for i in range(N):
        beta[t, i] = np.sum(A[i] * B[:, observations[t+1]] * beta[t+1])

# Step 3: Compute gamma (posterior probabilities)
P_obs = np.sum(alpha[T-1])
gamma = (alpha * beta) / P_obs

# Step 4: Compute xi (expected transitions)
xi = np.zeros((T-1, N, N))
for t in range(T-1):
    denom = 0
    for i in range(N):
        for j in range(N):
            denom += alpha[t, i] * A[i, j] * B[j, observations[t+1]] * beta[t+1, j]
    for i in range(N):
        for j in range(N):
            xi[t, i, j] = (alpha[t, i] * A[i, j] * B[j, observations[t+1]] * beta[t+1, j]) / denom

# Step 5: Update transition matrix A'
A_new = np.zeros((N, N))
for i in range(N):
    for j in range(N):
        A_new[i, j] = np.sum(xi[:, i, j]) / np.sum(gamma[:-1, i])

# Step 6: Update emission matrix B'
B_new = np.zeros((N, M))
for j in range(N):
    for k in range(M):
        mask = (observations == k)
        B_new[j, k] = np.sum(gamma[mask, j]) / np.sum(gamma[:, j])

# Step 7: Update initial probabilities
pi_new = gamma[0]

# Output all updated parameters
print("Updated Initial Probabilities (pi'):\n", pi_new)
print("\nUpdated Transition Matrix (A'):\n", A_new)
print("\nUpdated Emission Matrix (B'):\n", B_new)

print("\nPosterior Probabilities (gamma):\n", gamma)
print("\nExpected Transitions (xi):\n", xi)

To be continued in Part 2.