Two problems on writing a script to compute markov joint distribution (in python)

Question

I'm a new-learner of python, recently I'm working on some project to perform computation of Joint distribution of a markov process.

An example of a stochastic kernel is the one used in a recent study by Hamilton (2005), who investigates a nonlinear statistical model of the business cycle based on US unemployment data. As part of his calculation he estimates the kernel

pH :=   0.971 0.029 0
        0.145 0.778 0.077
        0     0.508 0.492

Here S = {x1, x2, x3} = {NG, MR, SR}, where NG corresponds to normal growth, MR to mild recession, and SR to severe recession. For example, the probability of transitioning from severe recession to mild recession in one period is 0.508. The length of the period is one month.

the excercise based on the above markov process is

With regards to Hamilton’s kernel pH, and using the same initial condition ψ = (0.2, 0.2, 0.6) , compute the probability that the economy starts and remains in recession through periods 0, 1, 2 (i.e., that xt = NG fort = 0, 1, 2).

My script is like

import numpy as np
## In this case, X should be a matrix rather than vector
## and we compute w.r.t P rather than merely its element [i][j]
path = []
def path_prob2 (p, psi , x2):  # X a sequence giving the path
     prob = psi                # initial distribution is an row vector
     for t in range(x2.shape[1] -1): # .shape[1] grasp # of columns
        prob = np.dot(prob , p)       # prob[t]: marginal distribution at period t
        ression = np.dot(prob, x2[:,t])
     path.append(ression)
     return path,prob



p = ((0.971, 0.029, 0    ),
      (0.145, 0.778, 0.077),
      (0    , 0.508, 0.492)) 
# p must to be a 2-D numpy array     
p = np.array(p)      

psi = (0.2, 0.2, 0.6)  
psi = np.array(psi)  
x2 = ((0,0,0),
      (1,1,1),
      (1,1,1))
x2 = np.array(x2)      
path_prob2(p,psi,x2)

During the execute process, two problems arise. The first one is , in the first round of loop, I don't need the initial distribution psi to postmultiply transaction matrix p, so the probability of "remaining in recession" should be 0.2+0.6 = 0.8, but I don't know how to write the if-statement. The second one is , as you may note, I use a list named path to collect the probility of "remaining in recession" in each period. And finally I need to multiply every element in the list one-by-one, I don't manage to find a method to implement such task , like path[0]*path[1]*path[2] (np.multiply can only take two arguments as far as I know). Please give me some clues if such method do exist. An additional ask is please give me any suggestion that you think can make the code more efficient. Thank you.

Answer 1

If I understood you correctly this should work (I'd love for some manual calculations for some of the steps/outcome), take notice of the fact that I didn't use if/else statement but instead started iterating from the second column:

import numpy as np

# In this case, X should be a matrix rather than vector
# and we compute w.r.t P rather than merely its element [i][j]
path = []


def path_prob2(p, psi, x2):  # X a sequence giving the path
    path.append(np.dot(psi, x2[:, 0]))  # first step
    prob = psi  # initial distribution is an row vector
    for t in range(1, x2.shape[1]):  # .shape[1] grasp # of columns
        prob = np.dot(prob, p)  # prob[t]: marginal distribution at period t
        path.append(np.prod(np.dot(prob, t)))
    return path, prob


# p must to be a 2-D numpy array
p = np.array([[0.971, 0.029, 0],
              [0.145, 0.778, 0.077],
              [0, 0.508, 0.492]])

psi = np.array([0.2, 0.2, 0.6])
x2 = np.array([[0, 0, 0],
               [1, 1, 1],
               [1, 1, 1]])

print path_prob2(p, psi, x2)

For your second question I believe Numpy.prod will give you a multiplication between all elements of a list/array.

You can use the prod as such:

>>> np.prod([15,20,31])
9300