### DEMDP02 
# Asset Replacement Model

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


```python
import numpy as np
import matplotlib.pyplot as plt
from compecon import BasisSpline, DPmodel, qnwnorm, demo, NLP
```

### 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$.



```python
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\}$. 


```python
dstates = [f'a={age+1:d}' for age in range(A)]
print('Discrete states are :', dstates)
```

    Discrete states are : ['a=1', 'a=2', 'a=3', 'a=4', 'a=5', 'a=6']


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


```python
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.


```python
options = ['keep', 'replace']
```

### 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$.



```python
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 
\begin{equation}p' = \bar{p} + \gamma(p-\bar{p}) + \epsilon \end{equation}
where $\epsilon \sim N(0,\sigma^2)$.


```python
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).


```python
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$.


```python
h = np.zeros((2, A),int)
h[0, :-1] = np.arange(1, A)
```

## Model Structure


```python
model = DPmodel(basis, profit, transition,
                i=dstates,
                j=options,
                discount=delta, e=e, w=w, h=h)
```

### SOLUTION
The ```solve``` method returns a pandas ```DataFrame```.


```python
S = model.solve()
S.head()
```

    Solving infinite-horizon model collocation equation by Newton's method
    iter change       time    
    ------------------------------
       0       3.7e+02    0.2599
       1       4.3e+01    0.5558
       2       1.2e+01    0.8437
       3       5.7e-01    1.1456
       4       3.7e-13    1.4555
    Elapsed Time =    1.46 Seconds





<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th></th>
      <th>unit profit</th>
      <th>i</th>
      <th>value</th>
      <th>resid</th>
      <th>j*</th>
      <th>value[keep]</th>
      <th>value[replace]</th>
    </tr>
    <tr>
      <th></th>
      <th>unit profit</th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th rowspan="5" valign="top">a=1</th>
      <th>0.000000</th>
      <td>0.000000</td>
      <td>0</td>
      <td>234.982583</td>
      <td>2.842171e-14</td>
      <td>keep</td>
      <td>234.982583</td>
      <td>206.124123</td>
    </tr>
    <tr>
      <th>0.001001</th>
      <td>0.001001</td>
      <td>0</td>
      <td>235.058576</td>
      <td>-6.159590e-07</td>
      <td>keep</td>
      <td>235.058576</td>
      <td>206.209036</td>
    </tr>
    <tr>
      <th>0.002001</th>
      <td>0.002001</td>
      <td>0</td>
      <td>235.134563</td>
      <td>-6.366804e-07</td>
      <td>keep</td>
      <td>235.134563</td>
      <td>206.293949</td>
    </tr>
    <tr>
      <th>0.003002</th>
      <td>0.003002</td>
      <td>0</td>
      <td>235.210544</td>
      <td>-2.317537e-07</td>
      <td>keep</td>
      <td>235.210544</td>
      <td>206.378862</td>
    </tr>
    <tr>
      <th>0.004002</th>
      <td>0.004002</td>
      <td>0</td>
      <td>235.286520</td>
      <td>4.292312e-07</td>
      <td>keep</td>
      <td>235.286520</td>
      <td>206.463774</td>
    </tr>
  </tbody>
</table>
</div>



## Analysis

### Plot Action-Contingent Value Functions


```python
# Compute and Plot Critical Unit Profit Contributions
pcrit = [NLP(lambda s: model.Value_j(s)[i].dot([1,-1])).broyden(0.0)[0] for i in range(A)]
vcrit = [model.Value(s)[i] for i, s in enumerate(pcrit)]
```


```python
fig1 = demo.figure('Action-Contingent Value Functions', 'Net Unit Profit','Value', figsize=[10,5])

cc = np.linspace(0.3,0.9,model.dims.ni)

for a, i in enumerate(dstates):
    plt.plot(S.loc[i,'value[keep]'] ,color=plt.cm.Blues(cc[a]), label='Keep ' + i)
    if pmin < pcrit[a] < pmax:
        demo.annotate(pcrit[a], vcrit[a], f'$p^*_{a+1}$', 'wo',(0, 0), fs=11, ms=18)
        print(f'Age {a+1:d}  Profit {pcrit[a]:5.2f}')
    
plt.plot(S.loc['a=1','value[replace]'], color=plt.cm.Oranges(0.5),label='Replace')
plt.legend()
```

### Plot Residual


```python
S['resid2'] = 100*S['resid'] / S['value']

fig2 = demo.figure('Bellman Equation Residual','Net Unit Profit','Percent Residual')
demo.qplot('unit profit','resid2','i',S)
plt.legend(model.labels.i)
```

### SIMULATION


```python
T = 50
nrep = 10000
sinit = np.full(nrep, pbar)
iinit = 0
data = model.simulate(T,sinit,iinit, seed=945)
data['age'] = data.i.cat.codes + 1
```

### Print Ergodic Moments


```python
frm = '\t{:<10s} = {:5.2f}'

print('\nErgodic Means')
print(frm.format('Price', data['unit profit'].mean()))
print(frm.format('Age', data['age'].mean()))
print('\nErgodic Standard Deviations')
print(frm.format('Price', data['unit profit'].std()))
print(frm.format('Age', data['age'].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


```python
data.head()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>time</th>
      <th>_rep</th>
      <th>i</th>
      <th>unit profit</th>
      <th>j*</th>
      <th>age</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>0</td>
      <td>0</td>
      <td>a=1</td>
      <td>1.0</td>
      <td>keep</td>
      <td>1</td>
    </tr>
    <tr>
      <th>1</th>
      <td>0</td>
      <td>1</td>
      <td>a=1</td>
      <td>1.0</td>
      <td>keep</td>
      <td>1</td>
    </tr>
    <tr>
      <th>2</th>
      <td>0</td>
      <td>2</td>
      <td>a=1</td>
      <td>1.0</td>
      <td>keep</td>
      <td>1</td>
    </tr>
    <tr>
      <th>3</th>
      <td>0</td>
      <td>3</td>
      <td>a=1</td>
      <td>1.0</td>
      <td>keep</td>
      <td>1</td>
    </tr>
    <tr>
      <th>4</th>
      <td>0</td>
      <td>4</td>
      <td>a=1</td>
      <td>1.0</td>
      <td>keep</td>
      <td>1</td>
    </tr>
  </tbody>
</table>
</div>




```python
subdata = data[data['_rep'].isin(range(3))]
```


```python
fig3 = demo.figure('Simulated and Expected Price', 'Period','Net Unit Profit')
demo.qplot('time','unit profit','_rep',subdata)
plt.plot(data[['time','unit profit']].groupby('time').mean(),'k--',label='mean')
plt.legend()
```

### Plot Expected Discrete State Path


```python
fig4 = demo.figure('Expected Machine Age', 'Period','Age')
demo.qplot('time','age','_rep',subdata)
plt.plot(data[['time','age']].groupby('time').mean(),'k--',label='mean')
plt.legend()
```

a8c023fd4d39193ab177687d6a130b81c899556f
ipynb
Jupyter Notebook
notebooks/dp/02 Asset Replacement Model.ipynb
daniel-schaefer/CompEcon-python
d3f66e04a7e02be648fc5a68065806ec7cc6ffd6
MIT
notebooks/dp/02 Asset Replacement Model.ipynb
daniel-schaefer/CompEcon-python
d3f66e04a7e02be648fc5a68065806ec7cc6ffd6
MIT
notebooks/dp/02 Asset Replacement Model.ipynb
daniel-schaefer/CompEcon-python
d3f66e04a7e02be648fc5a68065806ec7cc6ffd6
MIT
2021-06-01T03:47:35.000Z
2021-06-01T03:47:35.000Z
Qwen/Qwen-72B
1. YES
2. YES