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 \begin{equation} \pi_{t+1} = h(\pit, \epsilon{t+1}) \end{equation}
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 \begin{equation} j \in {0, 1} \end{equation}
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 \begin{equation} g(\pi, d, j, \epsilon) = h(\pi, \epsilon), j) \end{equation}
The reward function is \begin{equation} f(\pi, d, j) = \pi j − K_1(1 − d)j − K_0 d(1 − j) \end{equation}
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 \begin{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} \end{equation}
%matplotlib inline
from warnings import simplefilter
simplefilter('ignore')
from compecon import BasisSpline, DPmodel
from compecon.quad import qnwnorm
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 is \begin{equation} f(\pi, d, j) = \pi j − K_1(1 − d)j − K_0 d(1 − j) \end{equation}
where the exit cost is $K_0=0$ and the entry cost is $K_1=10$.
K0 = 0.0
K1 = 10
def profit(p, x, d, j):
return p * j - K1 * (1 - d) * j - K0 * d * (1 - j)
Assuming that the short-run profit $\pi$ is an exogenous Markov process
\begin{equation}
\pi_{t+1} = g(\pit,\epsilon{t+1}) = \bar{\pi} + \gamma(\pit − \bar{\pi}) + \epsilon{t+1}
\end{equation}
where $\bar{\pi}=1.0$ and $\gamma=0.7$.
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.
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
\begin{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\} \end{equation}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$.
n = 250
pmin = -20
pmax = 20
basis = BasisSpline(n, pmin, pmax, labels=['profit'])
print(basis)
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
\begin{equation}
h=\begin{bmatrix}0&0\1&1\end{bmatrix}
\end{equation}
model = DPmodel(basis,
profit, transition,
i=['idle', 'active'],
j=['idle', 'active'],
discount=0.9, e=e, w=w,
h=[[0, 0], [1, 1]])
To solve the model, we simply call the solve
method on the model
object.
S = mode
l.solve(print=True)
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()
We normalize the residuals as percentage of the value function. Notice the spikes at the "Profit entry" and "Profit exit" points.
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')
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.
T = 50
nrep = 50000
p0 = np.tile(pbar, (1, nrep))
d0 = 1
data = model.simulate(T, p0, d0, seed=945)
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()))
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()
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()