# Industry Entry-Exit Model¶

A firm operates in an uncertain profit environment. At the beginning of each period, the firm observes its potential short-run operating profit over the coming period $\pi$, which may be negative, and decides whether to operate, making a short-run profit $\pi$, or not operate, making a short-run profit $0$. Although the firm faces no fixed costs, it incurs a shutdown cost $K_0$ when it closes and a start-up cost $K_1$ when it reopens. The short-run profit $\pi$ is an exogenous continuous-valued Markov process \pi_{t+1} = h(\pit, \epsilon{t+1})

What is the optimal entry-exit policy? In particular, how low must the short-run profit be for an operating firm to close, and how high must the short-run profit be for nonoperating firm to reopen?

This is an infinite horizon, stochastic model with time $t$ measured in years. The state variables \begin{align} \pi &\in (−\infty,\infty)\ d &\in {0, 1} \end{align}

are the current short-run profit, a continuous variable, and the operational status of the firm, a binary variable that equals 1 if the firm is operating and 0 if the firm is not operating. The choice variable $$j \in {0, 1}$$

is the operating decision for the coming year, a binary variable that equals 1 if the firm operates and 0 if does not operate. The state transition function is $$g(\pi, d, j, \epsilon) = h(\pi, \epsilon), j)$$

The reward function is $$f(\pi, d, j) = \pi j − K_1(1 − d)j − K_0 d(1 − j)$$

The value of the firm, given that the current short-run profit is $\pi$ and the firm's operational status is $d$, satisfies the Bellman equation V(\pi,d) = \max_{j=0,1}\left{\pi j − K_1(1−d)j − K0d(1−j) + \delta E\epsilon V\left(h(\pi,\epsilon),j\right)\right}

In [1]:
%matplotlib inline
from warnings import simplefilter
simplefilter('ignore')
from compecon import BasisSpline, DPmodel
from demos.setup import np, plt, demo
import pandas as pd
from ggplot import ggplot, aes, geom_line, geom_point, theme_matplotlib


### The reward function¶

The reward function is $$f(\pi, d, j) = \pi j − K_1(1 − d)j − K_0 d(1 − j)$$

where the exit cost is $K_0=0$ and the entry cost is $K_1=10$.

In [2]:
K0     = 0.0
K1     = 10

def profit(p, x, d, j):
return p * j - K1 * (1 - d) * j - K0 * d * (1 - j)


### The transition function¶

Assuming that the short-run profit $\pi$ is an exogenous Markov process \pi_{t+1} = g(\pit,\epsilon{t+1}) = \bar{\pi} + \gamma(\pit − \bar{\pi}) + \epsilon{t+1}
where $\bar{\pi}=1.0$ and $\gamma=0.7$.

In [3]:
pbar   = 1.0
gamma  = 0.7

def transition(p, x, d, j, in_, e):
return pbar + gamma * (p - pbar) + e


In the transition function $\epsilon_t$ is an i.i.d. normal(0, $σ^2$), with $\sigma=1$. We discretize this distribution by using a discrete distribution, matching the first 10 moments of the normal distribution.

In [4]:
m = 5  # number of profit shocks
sigma  = 1.0
[e,w] = qnwnorm(m,0,sigma **2)


The collocation method calls for the analyst to select $n$ basis functions $\varphi_j$ and $n$ collocation nodes $(\pi_i,d_i)$, and form the value function approximant $V(\pi,d) \approx \sum_{j=1}^{n} c_j\varphi_j(\pi,d)$ whose coefficients $c_j$ solve the collocation equation

$$\sum_{j=1}^{n} c_j\varphi_j(\pi_i,d_i) = \max_{x\in\{0,1\}}\left\{\pi_i x − K_1(1−d_i)x − K_0 d_i(1−x) + \delta\sum_{k=1}^{m}\sum_{j=1}^{n}w_k c_j \varphi_j(\hat{\pi}_{ik},x)\right\}$$

where $\hat\pi_{ik}=g(\pi_i,\epsilon_k)$ and where $\epsilon_k$ and $w_k$ represent quadrature nodes and weights for the normal shock.

For the approximation, we use a cubic spline basis with $n=250$ nodes between $p_{\min}=-20$ and $p_\max=20$.

In [5]:
n = 250
pmin = -20
pmax =  20
basis = BasisSpline(n, pmin, pmax, labels=['profit'])
print(basis)

A 1-dimension Cubic spline basis:  using 250 Canonical nodes and 250 polynomials
___________________________________________________________________________
profit: 250 nodes in [-20.00,  20.00]

===========================================================================
WARNING! Class Basis is still work in progress


The Bellman equation is represeted by a DPmodel object, where we assume a discount factor of $\delta=0.9$. Notice that the discrete state transition is deterministic, with transition matrix $$h=\begin{bmatrix}0&0\1&1\end{bmatrix}$$

In [6]:
model = DPmodel(basis,
profit, transition,
i=['idle', 'active'],
j=['idle', 'active'],
discount=0.9, e=e, w=w,
h=[[0, 0], [1, 1]])


## SOLUTION¶

To solve the model, we simply call the solve method on the model object.

In [7]:
S = mode
l.solve(print=True)

Solving infinite-horizon model collocation equation by Newton's method
iter change       time
------------------------------
0       6.0e+01    0.0666
1       3.2e+01    0.1516
2       2.9e+00    0.2337
3       4.3e-01    0.3008
4       3.8e-14    0.3829
Elapsed Time =    0.38 Seconds


### Plot Action-Contingent Value Functions¶

In [8]:
S['ij'] = pd.Categorical(S.i.astype(np.ndarray) + '--' + S.j.astype(np.ndarray))
S.ij.cat.categories = ['Keep active firm open', 'Shut down', 'Reopen Firm', 'Remain idle']
fig=demo.qplot('profit', 'value_j', 'ij',
data=S[S.ij != 'Remain idle'],
geom='line',
xlab='Potential Profit',
ylab='Value of Firm')

#  Compute and Plot Critical Profit
ni, nj = model.dims['ni', 'nj']
vr = S['value_j'].reshape(ni, nj, -1)
sr = S['profit'].reshape(ni, nj, -1)[0, 0]
fig.draw()
offs = [(4, -6), (-8, 5)]

for i in range(2):
pcrit = np.interp(0, vr[i, 1] - vr[i, 0], sr)
vcrit  = np.interp(pcrit, sr, vr[i,0])
demo.annotate(pcrit, vcrit, '$p^*_{} = {:.2f}$'.format(i, pcrit), 'wo', offs[i], fs=11, ms=16,)
if i:
print('Profit Exit  = {:5.2f}'.format(pcrit))
else:
print('Profit Entry = {:5.2f}'.format(pcrit))
plt.show()

Profit Entry =  2.10
Profit Exit  = -2.30


### Plot Residual¶

We normalize the residuals as percentage of the value function. Notice the spikes at the "Profit entry" and "Profit exit" points.

In [9]:
S['resid2'] = 100 * (S.resid / S.value)
demo.qplot('profit','resid2','i',
data=S,
geom='line',
main='Bellman Equation Residual',
xlab='Potential Profit',
ylab='Percent Residual')

Out[9]:
<ggplot: (-9223371865891315789)>

## SIMULATION¶

We simulate the model 50000 times for a time horizon $T=50$, starting with an operating firm ($d=1$) at the long-term profit mean $\bar{\pi}$. To be able to reproduce these results, we set the random seed at an arbitrary value of 945.

In [10]:
T = 50
nrep = 50000
p0 = np.tile(pbar, (1, nrep))
d0 = 1
data = model.simulate(T, p0, d0, seed=945)

In [11]:
f = '\t{:21s} = {:5.2f}'
print('\nErgodic Means')
print(f.format('Profit Contribution', data['profit'].mean()))
print(f.format('Activity', (data['i'] == 'active').mean()))
print('\nErgodic Standard Deviations\n')
print(f.format('Profit Contribution', data['profit'].std()))
print(f.format('Activity', (data['i'] == 'active').std()))

Ergodic Means
Profit Contribution   =  1.00
Activity              =  0.94

Ergodic Standard Deviations

Profit Contribution   =  1.37
Activity              =  0.24


### Plot Simulated and Expected Continuous State Path¶

In [12]:
data2 = data[['time', 'profit']].groupby('time').mean()
data2['time'] = data2.index

ppp = ggplot(aes('time','profit','_rep'),
data=data[data._rep <3]) + \
geom_line() + theme_matplotlib()
ppp.draw()
plt.plot(data2.time,data2.profit)
plt.title('Simulated and Expected Profit Contribution')
plt.xlabel('Period')
plt.ylabel('Profit Contribution')
plt.show()


### Plot Expected Discrete State Path¶

In [13]:
subdata = data[['time', 'i']]
subdata['i'] = subdata['i'] == 'active'
subdata.groupby('time').mean().plot(legend=False)
plt.title('Probability of Operation')
plt.xlabel('Period')
plt.ylabel('Probability')

plt.show()