# Deterministic Optimal Economic Growth Model

**Randall Romero Aguilar, PhD**

This demo is based on the original Matlab demo accompanying the  <a href="https://mitpress.mit.edu/books/applied-computational-economics-and-finance">Computational Economics and Finance</a> 2001 textbook by Mario Miranda and Paul Fackler.

Original (Matlab) CompEcon file: **demdp06.m**

Running this file requires the Python version of CompEcon. This can be installed with pip by running

    !pip install compecon --upgrade

<i>Last updated: 2021-Oct-01</i>
<hr>

## About

Welfare maximizing social planner must decide how much society should consume and invest.  Model is of special interest because it has a known closed-form solution.

- States
    -     s       stock of wealth
- Actions
    -     k       capital investment
- Parameters
    -     beta    capital production elasticity
    -     delta   discount factor


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

### Model parameters

Assuming that the marginal productivity of capital is $\beta=0.5$ and the discount factor is $\delta=0.9$


```python
β, δ = 0.5, 0.9
```

## Analytic results

The steady-state values for this model are


```python
sstar = (β * δ) ** (β / (1 - β))   # steady-state wealth
kstar = β * δ * sstar                    # steady-state capital investment
vstar = np.log(sstar - kstar) / (1 - δ)     # steady-state value
pstar = 1 / (sstar * (1 - β * δ))        # steady-state shadow price
b = 1 / (1 - δ * β)

print('\n\nSteady-State')
for var, value in zip(['Wealth','Investment','Value','Shadow price'], [sstar,kstar,vstar,pstar]):
    print(f'\t{var:12s} = {value:8.4f}')
```

    
    
    Steady-State
    	Wealth       =   0.4500
    	Investment   =   0.2025
    	Value        = -13.9634
    	Shadow price =   4.0404


The true value function is
\begin{equation}
V(s) = v^* + \frac{1}{1-\delta\beta}\left(\log(s) -\log(s^*)\right)
\end{equation}


```python
def vtrue(wealth): # analytic value function
    return vstar + b * (np.log(wealth) - np.log(sstar))
```

The true policy function is
\begin{equation}
k(s) = \delta\beta s
\end{equation}


```python
def ktrue(wealth): #analytic policy function
    return δ*β*wealth
```

## Numeric results

### State space
The state variable is s="Wealth", which we restrict to $0\in[0.2, 1.0]$. 

Here, we represent it with a Chebyshev basis, with $n=15$ nodes.


```python
n, smin, smax = 15, 0.2, 1.0
basis = BasisChebyshev(n, smin, smax, labels=['Wealth'])
```

### Action space
The choice variable k="Investment" must be nonnegative.


```python
def bounds(s, i=None, j=None):
    return np.zeros_like(s), s[:]
```

### Reward function
The reward function is the utility of consumption=$s-k$.


```python
def reward(s, k, i=None, j=None):
    sk = s - k
    u = np.log(sk)
    ux= - sk ** -1
    uxx = - sk ** -2
    return u, ux, uxx
```

### State transition function
Next period, wealth will be equal to production from available initial capital $k$, that is $s' = k^\beta$


```python
def transition(s, k, i=None, j=None, in_=None, e=None):
    g = k ** β
    gx = β * k **(β - 1)
    gxx = (β - 1) * β * k ** (β - 2)
    return g, gx, gxx
```

### Model structure
The value of wealth $s$ satisfies the Bellman equation
\begin{equation*}
V(s) = \max_k\left\{\log(s-k) + \delta V(k^\beta)  \right\}
\end{equation*}

To solve and simulate this model,use the CompEcon class `DPmodel`


```python
growth_model = DPmodel(basis, reward, transition, bounds,
                       x=['Investment'],
                       discount=δ)
```

### Solving the model

Solving the growth model by collocation, using *Newton* algorithm and a maximum of 20 iterations


```python
options = dict(show=True,
               algorithm='newton',
               maxit=20)

snodes = growth_model.Value.nodes
S = growth_model.solve(vtrue(snodes), ktrue(snodes), **options)
```

    Solving infinite-horizon model collocation equation by Newton's method
    iter change       time    
    ------------------------------
       0       4.7e-07    0.0050
       1       1.6e-12    0.0080
    Elapsed Time =    0.01 Seconds


`DPmodel.solve` returns a pandas `DataFrame` with the following data:


```python
S.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>Wealth</th>
      <th>value</th>
      <th>resid</th>
      <th>Investment</th>
    </tr>
    <tr>
      <th>Wealth</th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0.200000</th>
      <td>0.200000</td>
      <td>-15.437865</td>
      <td>2.766030e-07</td>
      <td>0.090000</td>
    </tr>
    <tr>
      <th>0.205369</th>
      <td>0.205369</td>
      <td>-15.389699</td>
      <td>-2.099164e-07</td>
      <td>0.092416</td>
    </tr>
    <tr>
      <th>0.210738</th>
      <td>0.210738</td>
      <td>-15.342776</td>
      <td>-2.488242e-07</td>
      <td>0.094832</td>
    </tr>
    <tr>
      <th>0.216107</th>
      <td>0.216107</td>
      <td>-15.297033</td>
      <td>-1.102768e-07</td>
      <td>0.097248</td>
    </tr>
    <tr>
      <th>0.221477</th>
      <td>0.221477</td>
      <td>-15.252412</td>
      <td>5.644805e-08</td>
      <td>0.099664</td>
    </tr>
  </tbody>
</table>
</div>



We are also interested in the shadow price of wealth (the first derivative of the value function) and the approximation error.

To analyze the dynamics of the model, it also helps to compute the optimal change of wealth.


```python
S['shadow price'] = growth_model.Value(S['Wealth'],1)
S['error'] = S['value'] - vtrue(S['Wealth'])
S['D.Wealth'] = transition(S['Wealth'], S['Investment'])[0] - S['Wealth']
S.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>Wealth</th>
      <th>value</th>
      <th>resid</th>
      <th>Investment</th>
      <th>shadow price</th>
      <th>error</th>
      <th>D.Wealth</th>
    </tr>
    <tr>
      <th>Wealth</th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0.200000</th>
      <td>0.200000</td>
      <td>-15.437865</td>
      <td>2.766030e-07</td>
      <td>0.090000</td>
      <td>9.090758</td>
      <td>4.241867e-07</td>
      <td>0.100000</td>
    </tr>
    <tr>
      <th>0.205369</th>
      <td>0.205369</td>
      <td>-15.389699</td>
      <td>-2.099164e-07</td>
      <td>0.092416</td>
      <td>8.853206</td>
      <td>-3.099038e-08</td>
      <td>0.098631</td>
    </tr>
    <tr>
      <th>0.210738</th>
      <td>0.210738</td>
      <td>-15.342776</td>
      <td>-2.488242e-07</td>
      <td>0.094832</td>
      <td>8.627699</td>
      <td>-4.010325e-08</td>
      <td>0.097210</td>
    </tr>
    <tr>
      <th>0.216107</th>
      <td>0.216107</td>
      <td>-15.297033</td>
      <td>-1.102768e-07</td>
      <td>0.097248</td>
      <td>8.413362</td>
      <td>1.253504e-07</td>
      <td>0.095739</td>
    </tr>
    <tr>
      <th>0.221477</th>
      <td>0.221477</td>
      <td>-15.252412</td>
      <td>5.644805e-08</td>
      <td>0.099664</td>
      <td>8.209399</td>
      <td>3.151235e-07</td>
      <td>0.094220</td>
    </tr>
  </tbody>
</table>
</div>



### Solving the model by Linear-Quadratic Approximation

The `DPmodel.lqapprox` solves the linear-quadratic approximation, in this case arround the steady-state. It returns a LQmodel which works similar to the DPmodel object.

We also compute the shadow price and the approximation error to compare these results to the collocation results.


```python
growth_lq = growth_model.lqapprox(sstar, kstar)
L = growth_lq.solution(basis.nodes)
L['shadow price'] = L['value_Wealth']
L['error'] = L['value'] - vtrue(L['Wealth'])
L['D.Wealth'] = L['Wealth_next']- L['Wealth']
L.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>Wealth</th>
      <th>Investment</th>
      <th>value</th>
      <th>value_Wealth</th>
      <th>Wealth_next</th>
      <th>shadow price</th>
      <th>error</th>
      <th>D.Wealth</th>
    </tr>
    <tr>
      <th>Wealth</th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0.202191</th>
      <td>0.202191</td>
      <td>-0.020528</td>
      <td>-15.014820</td>
      <td>4.444953</td>
      <td>0.202192</td>
      <td>4.444953</td>
      <td>0.403234</td>
      <td>2.977733e-07</td>
    </tr>
    <tr>
      <th>0.219577</th>
      <td>0.219577</td>
      <td>-0.004880</td>
      <td>-14.937786</td>
      <td>4.416570</td>
      <td>0.219578</td>
      <td>4.416570</td>
      <td>0.330284</td>
      <td>2.769166e-07</td>
    </tr>
    <tr>
      <th>0.253590</th>
      <td>0.253590</td>
      <td>0.025731</td>
      <td>-14.788512</td>
      <td>4.361045</td>
      <td>0.253590</td>
      <td>4.361045</td>
      <td>0.217716</td>
      <td>2.361148e-07</td>
    </tr>
    <tr>
      <th>0.302742</th>
      <td>0.302742</td>
      <td>0.069968</td>
      <td>-14.576129</td>
      <td>4.280804</td>
      <td>0.302742</td>
      <td>4.280804</td>
      <td>0.107984</td>
      <td>1.771510e-07</td>
    </tr>
    <tr>
      <th>0.364886</th>
      <td>0.364886</td>
      <td>0.125897</td>
      <td>-14.313256</td>
      <td>4.179353</td>
      <td>0.364886</td>
      <td>4.179353</td>
      <td>0.031397</td>
      <td>1.026024e-07</td>
    </tr>
  </tbody>
</table>
</div>




```python
growth_lq.G
```




    array([[1.]])



## Plotting the results

### Optimal Policy


```python
fig1 = demo.figure('Optimal Investment Policy', 'Wealth', 'Investment')
plt.plot(S.Investment, label='Chebychev Collocation')
plt.plot(L.Investment, label='L-Q Approximation')
demo.annotate(sstar, kstar,'$s^*$ = %.2f\n$V^*$ = %.2f' % (sstar, kstar),'bo', (10, -17),ms=10)
plt.legend(loc= 'upper left')
```

### Value Function


```python
fig2 = demo.figure('Value Function', 'Wealth', 'Value')
plt.plot(S.value, label='Chebychev Collocation')
plt.plot(L.value, label='L-Q Approximation')

demo.annotate(sstar, vstar, f'$s^* = {sstar:.2f}$\n$V^* = {vstar:.2f}$', 'bo', (10, -17),ms=10)
plt.legend(loc= 'upper left')
```

### Shadow Price Function


```python
fig3 = demo.figure('Shadow Price Function', 'Wealth', 'Shadow Price')
plt.plot(S['shadow price'], label='Chebychev Collocation')
plt.plot(L['shadow price'], label='L-Q Approximation')
demo.annotate(sstar, pstar,f'$s^* = {sstar:.2f}$\n$\lambda^* = {pstar:.2f}$', 'bo', (10, 17),ms=10)
plt.legend(loc= 'upper right')
```

### Chebychev Collocation Residual and Approximation Error vs. Linear-Quadratic Approximation Error


```python
fig4 = plt.figure(figsize=[12, 6])
demo.subplot(1, 2, 1, 'Chebychev Collocation Residual\nand Approximation Error', 'Wealth', 'Residual/Error')
plt.hlines(0,smin,smax,'k',linestyles='--')
plt.plot(S[['resid', 'error']])
plt.legend(['Residual','Error'], loc='lower right')
plt.ticklabel_format(style='sci', axis='y', scilimits=(-1,1))

demo.subplot(1, 2, 2, 'Linear-Quadratic Approximation Error', 'Wealth', 'Error')
plt.hlines(0,smin,smax,'k',linestyles='--')
plt.plot(L['error'], label='Error')
plt.legend(loc='upper left')
plt.ticklabel_format(style='sci', axis='y', scilimits=(-1,1))
```

### Wealth dynamics

Notice how the steady-state is stable in the Chebyshev collocation solution, but unstable in the linear-quadratic approximation. In particular, simulated paths of wealth in the L-Q approximation will converge to zero, unless the initial states is within a small neighborhood of the steady-state.


```python
fig5 = demo.figure('Wealth dynamics', 'Wealth', 'Wealth change', figsize=[8,5])
plt.plot(S['D.Wealth'], label='Chebychev Collocation')
plt.plot(L['D.Wealth'], label='L-Q Approximation')
plt.hlines(0,smin,smax,linestyles=':')

demo.annotate(sstar, 0, f'$s^* = {sstar:.2f}$\n$\Delta s^* = {0:.2f}$', 'bo', (10, 10),ms=10,fs=11)
plt.legend(loc= 'lower left')
```

## Simulating the model

We simulate 20 periods of the model starting from $s=s_{\min}$


```python
T = 20
data = growth_model.simulate(T, smin)
```

### Simulated State and Policy Paths


```python
opts = dict(spec='r*', offset=(-5, -5), fs=11, ha='right')

fig6 = demo.figure('State and Policy Paths','Period', 'Wealth/Investment',[0, T + 0.5])
plt.plot(data[['Wealth', 'Investment']])
demo.annotate(T, sstar, 'steady-state wealth\n = %.2f' % sstar, **opts)
demo.annotate(T, kstar, 'steady-state investment\n = %.2f' % kstar, **opts)
```


```python
#demo.savefig([fig1,fig2,fig3,fig4,fig5,fig6])
```

ad35044efdffaad60fe2f70d13b41e68101e5051
ipynb
Jupyter Notebook
_build/jupyter_execute/notebooks/dp/06 Deterministic Optimal Economic Growth Model.ipynb
randall-romero/CompEcon-python
c7a75f57f8472c972fddcace8ff7b86fee049d29
MIT
2016-12-14T13:21:27.000Z
2020-08-23T21:04:34.000Z
_build/jupyter_execute/notebooks/dp/06 Deterministic Optimal Economic Growth Model.ipynb
randall-romero/CompEcon
c7a75f57f8472c972fddcace8ff7b86fee049d29
MIT
2017-09-10T04:48:54.000Z
2018-03-31T01:36:46.000Z
_build/jupyter_execute/notebooks/dp/06 Deterministic Optimal Economic Growth Model.ipynb
randall-romero/CompEcon-python
c7a75f57f8472c972fddcace8ff7b86fee049d29
MIT
2017-02-25T08:10:38.000Z
2020-05-15T09:49:16.000Z
Qwen/Qwen-72B
1. YES
2. YES