# Asset Replacement Model¶

Profit-maximizing entrepreneur must decide when to replace an aging asset.

In [1]:
%matplotlib inline
from warnings import simplefilter
simplefilter('ignore')
from compecon import BasisSpline, DPmodel
from demos.setup import np, plt, demo


### Model Parameters¶

The maximum age of the asset is $A=6$, production function coefficients are $\alpha =[50.0, -2.5, -2.5]$, net replacement cost is $\kappa = 40$, long-run mean unit profit is $\bar{p} = 1$, the unit profit autoregression coefficient is $\gamma = 0.5$, the standard deviation of unit profit shock is $\sigma = 0.15$, and the discount factor is $\delta = 0.9$.

In [2]:
A       = 6
alpha   = np.array([50, -2.5, -2.5])
kappa   = 40
pbar    = 1.0
gamma   = 0.5
sigma   = 0.15
delta   = 0.9


### State space¶

The state variables are the current unit profit $p$ (continuous) and the age of asset $a\in\{1,2,\dots,A\}$.

Here, we approximate the unit profit with a cubic spline basis of $n=200$ nodes, over the interval $p\in[0,2]$.

In [3]:
n  = 200
pmin, pmax = 0.0, 2.0
basis = BasisSpline(n, pmin, pmax, labels=['unit profit'])


### Action space¶

There is only one choice variable $j$: whether to keep(0) or replace(1) the asset.

### Reward Function¶

The instant profit depends on the asset age and whether it is replaced. An asset of age $a=A$ must be replaced. If the asset is replaced, the profit is $50p$ minus the cost of replacement $\kappa$; otherwise the profit depends on the age of the asset, $(\alpha_0 + \alpha_1 a + \alpha_2 a^2)p$.

In [4]:
def profit(p, x, i, j):
a = i + 1
if j or a == A:
return p * 50 - kappa
else:
return p * (alpha[0] + alpha[1] * a + alpha[2] * a ** 2 )


### State Transition Function¶

The unit profit $p$ follows a Markov Process $$p' = \bar{p} + \gamma(p-\bar{p}) + \epsilon$$ where $\epsilon \sim N(0,\sigma^2)$.

In [5]:
def transition(p, x, i, j, in_, e):
return pbar + gamma * (p - pbar) + e


The continuous shock must be discretized. Here we use Gauss-Legendre quadrature to obtain nodes and weights defining a discrete distribution that matches the first 10 moments of the Normal distribution (this is achieved with $m=5$ nodes and weights).

In [6]:
m = 5
e, w = qnwnorm(m,0,sigma ** 2)


On the other hand, the age of the asset is deterministic: if it is replaced now it will be new ($a=0$) next period; otherwise its age will be $a+1$ next period if current age is $a$.

In [7]:
h = np.zeros((2, A))
h[0, :-1] = np.arange(1, A)


## Model Structure¶

In [8]:
model = DPmodel(basis, profit, transition,
i=[a + 1 for a in range(A)],
j=['keep', 'replace'],
discount=delta, e=e, w=w, h=h)


### SOLUTION¶

The solve method returns a pandas DataFrame.

In [9]:
S = model.solve()

Solving infinite-horizon model collocation equation by Newton's method
iter change       time
------------------------------
0       3.7e+02    0.2287
1       4.3e+01    0.5015
2       1.2e+01    0.7568
3       5.7e-01    0.9986
4       4.9e-13    1.2563
Elapsed Time =    1.26 Seconds


## Analysis¶

### Plot Action-Contingent Value Functions¶

In [10]:
fig = demo.qplot('unit profit', 'value_j', 'i',
data=S,
main='Action-Contingent Value Functions',
xlab='Net Unit Profit',
ylab='Value')

ni, nj = model.dims['ni', 'nj']
vr = S['value_j'].reshape(ni, nj, -1)
pr = S['unit profit'].reshape(ni, nj, -1)[0, 0]
fig.draw()

print('Critical Replacement Profit\n')
for a in range(A-1):
pcrit = np.interp(0.0, vr[a, 1] - vr[a, 0], pr, np.nan, np.nan)
vcrit = np.interp(pcrit, pr, vr[a, 0])
if np.isnan(pcrit):
continue

demo.annotate(pcrit, vcrit, '$p^*_' + str(a+1) + '$', 'wo',
(0, 0), fs=11, ms=18)
print('   Age {:2d}  Profit {:5.2f}'.format(a, pcrit))

plt.xlim(pmin,pmax)
plt.show()

Critical Replacement Profit

Age  1  Profit  1.50
Age  2  Profit  0.66
Age  3  Profit  0.38
Age  4  Profit  0.25


### Plot Residual¶

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

Out[11]:
<ggplot: (163843090040)>

### SIMULATION¶

In [12]:
T = 50
nrep = 10000
sinit = np.full(nrep, pbar)
iinit = 0
data = model.simulate(T,sinit,iinit, seed=945)

In [13]:
frm = '\t{:<10s} = {:5.2f}'

print('\nErgodic Means')
print(frm.format('Price', data['unit profit'].mean()))
print(frm.format('Age', data.i.mean()))
print('\nErgodic Standard Deviations')
print(frm.format('Price', data['unit profit'].std()))
print(frm.format('Age', data.i.std()))

Ergodic Means
Price      =  1.00
Age        =  2.01

Ergodic Standard Deviations
Price      =  0.17
Age        =  0.83


### Plot Simulated and Expected Continuous State Path¶

In [14]:
demo.qplot('time', 'unit profit', '_rep',
data=data[data['_rep'] < 3],
geom='line',
main='Simulated and Expected Price',
ylab='Net Unit Profit',
xlab='Period')

Out[14]:
<ggplot: (-9223371873012094839)>

### Plot Expected Discrete State Path¶

In [15]:
data[['time', 'i']].groupby('time').mean().plot(legend=False)
plt.title('Expected Machine Age')
plt.xlabel('Period')
plt.ylabel('Age')

Out[15]:
<matplotlib.text.Text at 0x2625e196828>