[Economic Theory]Hopenhayn Model: Python

In continuation with the last post, in this one, we will be providing a step by step handbook to program and solve for the model that we looked at.

If you haven’t read it yet, please read the following first.

[Economic Theory] Hopenhayn Model

Import packages

import numpy as np
from scipy.stats import norm 
from scipy.stats import lognorm
import matplotlib.pyplot as plt

from math import sqrt
from math import exp
from math import log

from statistics import stdev
from sklearn.preprocessing import StandardScaler

Define Parameters

1. For AR(1) process

rho=0.9415
sigma=0.1908

2. Entrant profitability distribution

mu_e=-0.44 

3. Upper bar for fixed cost

phi_f_up=2.3

4. Firm’s optimization

alpha=2/3
beta=1/1.04 #4 percent interest rate

5. Equilibrium conditions

Free entry condition

Market clearing

N_bar=0.597 #employment rate
phi_e=5.18 

C&P)

rho=0.9415
sigma=0.1908
mu_e=-0.44
phi_f_up=2.3
alpha=2/3
beta=1/1.04 
N_bar=0.597 
phi_e=5.18

Discretize Random variable log z

Gaussian AR(1) process

There are many ways to discritize time series random variable. Tauchen is just one way to do it.

1. Make a discrete equi-spaced grid for log z.

Under the covariance stationarity condition of AR(1) process, the unconditional standard deviation of log(z) can be derived.

By the stationarity, we can rewrite as following.

Then, we can define a maximum (minimum) of the grid by the multiple m>0 (-m) of the unconditional standard deviation. Then, construct a equi-spaced grid with s number of points. In our project, we let m=4, s=101.

def log_z_grid(s):
  m=4
  std=sqrt(sigma**2/(1-rho**2))
  grid=np.linspace(-m*std,m*std,s)
  return grid

s=101
log_z=log_z_grid(s)

2. Construct a transition matrix (markov chain) F[z’|z] of size s by s.

For points from index 1 to s, first let’s consider the ones in the middle j=2,…,s-1. Starting from any point z_i, our goal is to get a probability that it goes to a point z_j. i and j corresponds to row number and column number respectively. We define that a point z_j is reached if it arrives within the distance less than half of the width between points.

Then for the points at the very end, since there is only one way to go,

Also, since we not only need the probability itself (going to one z_j) but also the cumulative probability (going to any z below z_j) we define them f_markov and F_markov, respectively.

def F_markov_cum_dis(f_markov):
  F_markov=np.zeros(f_markov.shape)
  F_markov[:,0]=f_markov[:,0]
  for i in range(1,len(f_markov)):
    F_markov[:,i]=f_markov[:,i]+F_markov[:,i-1]
  return F_markov

def tauchen(log_z_array):
  w=log_z_array[1]-log_z_array[0]
  s=len(log_z_array)
  f_markov=np.zeros((s,s))
  for i in range(s): #s=len(z_grid)
    for j in range(s):
      if j == np.min(range(s)): #from any i arrives first in grid
        f_markov[i,j]=norm.cdf((log_z_array[j]-rho*log_z_array[i]+w/2)/sigma)
      elif j == np.max(range(s)): #from any i arrives last in grid 
        f_markov[i,j]=1-norm.cdf((log_z_array[j]-rho*log_z_array[i]-w/2)/sigma)
      else:
        f_markov[i,j]=norm.cdf((log_z_array[j]-rho*log_z_array[i]+w/2)/sigma)-\
                        norm.cdf((log_z_array[j]-rho*log_z_array[i]-w/2)/sigma)    
  F_markov=F_markov_cum_dis(f_markov)  
  return f_markov, F_markov

We can plot some rows of markov chain (f_markov) and cumulative distribution (F_markov), e.g., starting from z_0, z_40, etc.

Transition matrix

Cumulative

Gaussian distribution for entrants

Using the same grid points for log z, since we assume that F_E(z) is normal distribution, it is easy to get both pdf and cdf.

def cumulative(pdf):
  cdf=np.zeros((len(pdf),1))
  cdf[0]=pdf[0]
  for i in range(1,len(pdf)):
    cdf[i]=cdf[i-1]+pdf[i]
  return cdf #column vector

def entrant(log_z_array):
  w=log_z_array[1]-log_z_array[0]
  pdf=[]
  for i in range(len(log_z_array)):
    prob=norm.cdf((log_z_array[i]+w/2),mu_e,sigma)-norm.cdf((log_z_array[i]-w/2),mu_e,sigma)
    pdf.append(prob)
  pdf=np.array(pdf)
  return pdf

f_en=entrant(log_z)
F_en=cumulative(f_en)

Algorithm Roadmap

Now that all the ingredients are ready (parameters and discretized profitability), we will briefly explain how the algorithm would look like.

1. Outer loop on equilibrium objects: given wage and entry mass (M_E)
2. Inner loop 1. Firm’s optimization problem
   Given wage find a value function that solves for the Bellman equation. Do it until convergence of value function.
3. Inner loop 2. Stationarity on firm distribution and mass
   Given entry mass, guess initial operating firm distribution and mass of operating firms. Then, solve for the two stationarity conditions under the mass and distribution of operating firms converge.
4. Convergence on outer loop: free entry condition and labor market clearing condition

Here we can separate the process by

1. Guess for wage, solve for value function and compute value upon entry (V_E(z)) then check for free entry condition. Update the guess on wage until free entry condition is satisfied.
2. Guess for entry mass then solve for two conditions on firm distribution and mass.
3. Check for market clearing condition.

By separating the process, we can be free of making guesses of two unknowns, rather once equilibrium wage is found we are sure that it is correct one.

Value Function Iteration

We can rewrite the bellman equation that we saw in the last post in terms of iterations k. Under certain conditions (which hold in this case and in most of macro models), the bellman is contraction mapping. In other words, starting from any guess on value function denoted as $V^{0}(z)$ and usually a zero vector of size (sx1), we are guaranteed to converge to a true value function.

We continue the iteration over value function until

Static choice of n

As we discussed in the previous post, we can separate the optimal choices of two choice variables: labour demand and exit threshold. By FOC result given z and W, optimal choice of labour demand is

z=np.exp(log_z)

def profit_func(wage,z):
  z=z.reshape((len(z),1))
  best_n_array=(wage/(alpha*z))**(1/(alpha-1))
  max_profit=np.multiply(z,best_n_array**alpha)-wage*best_n_array
  max_profit=max_profit.reshape((len(max_profit),))
  return best_n_array, max_profit

Optimal exit threshold phi_star.

In the previous post, we discussed that the threshold value for the fixed cost is that makes the firm in different between exit (zero profit) and discounted expected value that it gains if move to next period.

def phi_star_func(value,f_markov):
  s=len(value)
  value=value.ravel()
  f=np.multiply(value,f_markov)
  phi_star=beta*(1-delta)*np.sum(f,axis=1)
  return phi_star

Also, given this phi_star, we can compute the expected fixed cost that this firm has to pay next period.

def phi_exp_func(phi_star):
  s=len(phi_star)
  phi_exp=np.zeros((s,))
  for i in range(s):
    if phi_star[i]<=phi_f_up:
      phi_exp[i]=((phi_star[i]**2)/(2*phi_f_up))
    else:
      phi_exp[i]=phi_f_up/2
  return phi_exp

Here note that there is upper bar in the uniform distribution so that we have to differentiate the cases.

def phi_exp_func(phi_star):
  s=len(phi_star)
  phi_exp=np.zeros((s,))
  for i in range(s):
    if phi_star[i]<=phi_f_up:
      phi_exp[i]=((phi_star[i]**2)/(2*phi_f_up))
    else:
      phi_exp[i]=phi_f_up/2
  return phi_exp

VFI and Free entry condition: Guess wage

Now, we have constructed 3 parts required to update the value function. Given W and value function V^{(k)}, we can compute: max_profit+ phi_star- phi_exp which will become updated value function for iteration k+1.

def VFI(max_profit,z,f_markov,verbose=False):
  max_iter=3000
  tol_VFI=1e-10
  #To keep track of iterations. 
  v=np.zeros((len(z),max_iter)) 
  phi_star=np.zeros((len(z),max_iter))

  v_0=np.zeros((len(z),))#Initial guess
  phi_star_0= phi_star_func(value=v_0,f_markov=f_markov)
  phi_star[:,0]=phi_star_0.ravel()

  for i in range(1,max_iter):
    v[:,i]=max_profit+phi_star[:,i-1]-phi_exp_func(phi_star[:,i-1])
    phi_star[:,i]=phi_star_func(value=v[:,i],f_markov=f_markov)
    if np.max(abs(v[:,i]-v[:,i-1]))<tol_VFI:
      break

  if verbose:
    print('VFI iteration=',i)
  return v[:,i], phi_star[:,i], v[:,:i+1], phi_star[:,:i+1]

Then given the converged value function, we can compute value upon entry.

#given the result of VFI, compute value upon entry
def V_E_func(value_conv,f_en): 
  value_conv=value_conv.ravel()
  f_en=f_en.ravel()
  V_E=np.sum(np.multiply(value_conv,f_en))
  return V_E

How to use these functions:

opt_wage=0.399
best_n, max_profit=profit_func(wage=opt_wage,z=z)
value_conv,phi_conv,v_total,phi_total=VFI(max_profit,z=z,f_markov=f_markov)
V_E=V_E_func(value_conv=value_conv,f_en=f_en)
print(abs(V_E-phi_e))

Note that we can construct a function that inputs wage and outputs the absolute difference between V_E and phi_e.

By using, zero finder of the function, we were able to get optimal wage at 0.399.

Optimal wage

To plot:

fig, axs = plt.subplots(2,figsize=(6,6),sharex=True)
fig.suptitle('wage={:.3f}'.format(opt_wage))
axs[0].plot(z,profit_func(opt_wage,z=z)[0])
axs[1].plot(z,profit_func(opt_wage,z=z)[1])
axs[0].set_title("best_n")
axs[1].set_title("profit")
axs[1].set(xlabel="z")
plt.show()

def VFI_check(v_total,phi_f_total):
  fig, axs = plt.subplots(2,figsize=(6,8),sharex=True)
  # fig.suptitle('Convergence check')
  iter=list(range(0,v_total.shape[1],5))
  iter.append(v_total.shape[1]-1)
  for k in iter:
    axs[0].plot(z,v_total[:,k])
    axs[1].plot(z,phi_f_total[:,k])
    axs[0].set_title("V(z)")
    axs[1].set_title("$\phi^{*}(z)$")
    axs[1].set(xlabel="z")
  return

Firm Distribution

As we discussed earlier, we will make guesses on entrant firm mass operating firm mass and distribution. We can rewrite the two equations in terms of iterations.

Mass of firm

• $G(\phi^{*}_{F}(z))$: given z, probability of drawing a fixed cost so that the firm can operate next period.
• $f_o^{(k)}(z)$: probability that a firm has shock z. (weights or probabilities)

→ The weighted sum would be expected probability of any firm to stay in the market.

def G_phi_star_func(phi_conv):
  s=len(phi_conv)
  G_phi_star=np.zeros((s,))
  for i in range(s):
    if phi_conv[i]<=phi_f_up:
      G_phi_star[i]=phi_conv[i]*(1/phi_f_up)
    else:
      G_phi_star[i]=1
  return G_phi_star

Distribution of firm

• $G(\phi_F^{*}(z))$: Given z, probability to draw fixed cost below threshold → ith element
• $F[z’|z]$: Given z, probability to go any shock below z’ tmw →ith row
• $f_o^{(k)}(z)$: probability to be in z → ith element

→ Expected probability that a firm goes to next period and faces shock less than z’.

We will keep updating the two equations until

def firm_dist(M_op_0, M_en_0,F_markov,phi_conv,log_z,verbose=False):
  G_phi_star=G_phi_star_func(phi_conv)
  F_en=cumulative(entrant(log_z)).ravel()
  max_iter=10
  tol_F=1.e-10
  s=len(phi_conv)

  M_op=np.zeros((1,max_iter))
  f_op=np.zeros((s,max_iter))
  F_op=np.zeros((s,max_iter))
  #initial guesses
  M_op[:,0]=M_op_0
  f_op[:,0]=entrant(log_z) #same as entrants. 
  F_op[:,0]=cumulative(f_op[:,0]).ravel()


  arg=np.multiply(G_phi_star.reshape((s,1)),F_markov) #instead of nested np.multiply define as arg

  for it in range(1,max_iter):
    M_op[:,it]=(1-delta)*M_op[:,it-1]* np.sum(np.multiply(G_phi_star,f_op[:,it-1].ravel()))+M_en_0
    matrix_f=arg*f_op[:,it-1].reshape((s,1))
    integral=np.sum(matrix_f,axis=0)
    F_op[:,it]=(1-delta)*(M_op[:,it-1]/M_op[:,it])*integral +(M_en_0/M_op[:,it])*(F_en)
    f_op[:,it]=decumulative(F_op[:,it]).ravel()

    if ((np.max(abs(F_op[:,it]-F_op[:,it-1])) <tol_F) and (abs(M_op[:,it]-M_op[:,it-1])< tol_F)):
      break
  if verbose:
    print('Firm iterations=',it)
  return M_op[:,it], F_op[:,it], f_op[:,it], M_op[:,:it+2], F_op[:,:it+2], f_op[:,:it+2]

Labour Market Clearing

To be sure that the model solved for the equilibrium we have to check for the labour market clearing condition. Note that this requires both the result from VFI (more precisely optimal labour demand from guess of wage) and converged mass and distribution of operating firm.

def employment_rate(best_n,M_op_con,f_op_con):
  N=M_op_con*np.sum(np.multiply(best_n.ravel(),f_op_con))
  return N[0]

Unlike the free entry condition, this requires guess on both M_E and M_O, we computed employment rate and compared it with N_bar (which we were able to get precisely). We computed the distance from the N_bar given different combinations of guesses and got the one that was best.

Note here that unlike VFI, this procedure is not a contraction mapping so that there are multiple possible routes that the iterations could lead. In that sense, on top of employment rate, we wanted to have similar average exit rate and employees per firm as data given in the paper so we also included that in terms of choosing the initial guesses.

def exit_rate(phi_conv,f_op_con):
  G_phi_star=G_phi_star_func(phi_conv=phi_conv)
  survival_rate=(1-delta)*np.sum(np.multiply(G_phi_star,f_op_con.ravel()))
  return (1-survival_rate)

def firm_employment(best_n,f_op_con):
  firm_level_emp=np.sum(np.multiply(best_n.ravel(),f_op_con))
  return firm_level_emp

for M_en_0 in M_en_0_list:
  for M_op_0 in M_op_0_list:
    M_op_con,F_op_con,f_op_con,_,_,_=firm_dist(M_op_0=M_op_0, M_en_0=M_en_0,F_markov=F_markov,log_z=log_z,phi_conv=phi_conv)
    N=employment_rate(best_n,M_op_con,f_op_con)
    # print(M_op_con,N)
    if abs(N-0.597)<0.001:
      firm_N=firm_employment(best_n,f_op_con)
      e=exit_rate(phi_conv=phi_conv,f_op_con=f_op_con)
      print('when M_en_0 and M_op_0=',[M_en_0,M_op_0],'N,firm_N,e',N,firm_N,e,'M_op_con',M_op_con)

To check for convergence:

def FD_check_z(f_total,F_total):
  fig, axs = plt.subplots(2,figsize=(6, 7))
  # fig.suptitle('Convergence check')
  for k in range(f_total.shape[1]):
    axs[0].plot(log_z,f_total[:,k])
    axs[1].plot(log_z,F_total[:,k])
    axs[0].set_title("Density")
    axs[1].set_title("Distribution")
    axs[1].set(xlabel="log(z)")
    axs[0].annotate("",xy=(0.5,0.02),xytext=(0,0.08),fontsize=15, arrowprops=dict(
                             color='red',
                             ls='--'))
    axs[1].annotate("",xy=(0,0.6),xytext=(-0.3,0.8),fontsize=15, arrowprops=dict(
                             color='red',
                             ls='--'))
  plt.show()
  return

Thank you for reading. Hope that this post can help you how to approach solving different macro models!

[Citations]

Tauchen, G. (1986). Finite State Markov-chain Approximations to Univariate and Vector Autoregressions. Economics Letters 20 (2), 177–181.