<jupyter_start><jupyter_text>Your First Consumption-Saving Model Consider the following **consumption-saving model**:$$ \begin{aligned}v_{t}(m_{t}) &= \max_{c_t}\frac{c^{1-\rho}}{1-\rho} + \beta \mathbb{E}_t[v_{t+1}(m_{t+1})] \\& \text{s.t.} \\ a_{t}&=m_{t}-c_{t} \\m_{t+1}&= R a_{t} + \tilde{\xi}_{t+1}\\a_t&\geq 0\\\end{aligned}$$where$$ \begin{aligned}\log\xi_{t+1} \sim \mathcal{N}(-0.5\sigma_{\xi}^{2},\sigma_{\xi}^{2})\end{aligned} $$In the **last period** everything is consumed such that$$ \begin{aligned}v_{T}(m_{T}) = \frac{m_T^{1-\rho}}{1-\rho} \end{aligned} $$ To **solve this model using the ConSav package** we need to define the following functions:1. `settings(model)` Choose fundamental settings2. `setup(model)`: Set baseline parameters3. `allocate(model)`: Create grids and allocate solution and simulation memory4. `solve(model)`: Solve the model5. `simulate(model)`: Simulate the modelThen we collect all of these functions as methods in a **class**.  Numerical integration **General problem:** How can we calculate $$\mathbb{E}(f(x))=\int f(x)g(x)dx$$where 1. $f:\mathbb{R}\rightarrow\mathbb{R}$ is some function2. $g(x)$ is the probability distribution function (PDF) for $x$ and $G(x)$ is the CDF**General solution:** Turn it into a discrete sum$$\mathbb{E}(f(x))\approx\sum_{i=1}^{S}\omega_{i}f(x_{i})$$**Gauss-Hermite:** Specific clever way of choosing the $\omega_i$'s and $x_i$'s. Provided for a log-normal variable by ConSav in `normal_gauss_hermite`.  Setup<jupyter_code>%matplotlib inline
%load_ext autoreload
%autoreload 2

import time
import numpy as np
from numba import njit, prange
from scipy import optimize

import matplotlib.pyplot as plt
plt.style.use('seaborn-whitegrid')
prop_cycle = plt.rcParams['axes.prop_cycle']
colors = prop_cycle.by_key()['color']

# consav, install with pip install ConSav
from consav import ModelClass, jit, linear_interp, golden_section_search
from consav.misc import elapsed
from consav.grids import nonlinspace
from consav.quadrature import log_normal_gauss_hermite<jupyter_output><empty_output><jupyter_text>Components The **setings** function:1. `namespaces`: Will be availible as `.par`, `.sol` and `.sol` and can be used in numba jitted functions.2. `other_attrs`: Other attributes to be saved (cannot be used in numba jitted functions).3. `savefolder`: Where to save the model.4. `not_floats`: List of non-float elements in namespaces (ensures safe type inference).<jupyter_code>def settings(model):
    """ choose fundamental settings """
    
    # a. namespaces
    model.namespaces = ['par','sol','sim']
    
    # b. other attributes
    model.other_attrs = []
    
    # c. savefolder
    model.savefolder = 'saved'
    
    # b. list not-floates to safe type inferences
    model.not_floats = ['T','Nxi','Nm','simN'] 
    # rho = 2 e.g. raises an error, while T = 2 does not<jupyter_output><empty_output><jupyter_text>The **setup** function:<jupyter_code>def setup(model):
    """ set baseline parameters (independent elements) """
        
    par = model.par

    # b. parameters
    par.T = 100 # number of periods
    par.rho = 2.0 # CRRA parameter
    par.beta = 0.96 # discount factor
    par.sigma_xi = 0.1 # std. of transitory shock
    par.Nxi  = 8 # number of quadrature points for xi
    par.R = 1.04 # return factor
    par.Nm = 100 # number of nodes in grid for m
    par.m_min = 1e-6 # minimum point in grid for m
    par.m_max = 20.0 # maximum point in grid for m
    par.m_phi = 1.3 # curvature parameter for grid for m
    par.simN = 10_000 # number of persons in simulation<jupyter_output><empty_output><jupyter_text>The **allocate** function:<jupyter_code>def allocate(model):
    """ allocate model, i.e. create grids and allocate solution and simluation arrays (dependent elements) """
    
    par = model.par
    sol = model.sol
    sim = model.sim   
    
    # a. grids
    
    # non-equally spaced grid
    par.grid_m = nonlinspace(par.m_min,par.m_max,par.Nm,par.m_phi)
    
    # nodes and weights using Gauss-Hermite
    par.xi, par.xi_w = log_normal_gauss_hermite(sigma=par.sigma_xi,n=par.Nxi) 
    
    # b. solution
    sol.c = np.zeros((par.T,par.Nm))
    sol.v = np.zeros((par.T,par.Nm))
    
    # c. simulation
    sim.m = np.zeros((par.T,par.simN))
    sim.c = np.zeros((par.T,par.simN))
    sim.a = np.zeros((par.T,par.simN))
    sim.xi = np.zeros((par.T,par.simN))<jupyter_output><empty_output><jupyter_text>The **solve** function:<jupyter_code>def value_of_choice(c,t,m,model):
    """ value of choice """
    
    par = model.par
    sol = model.sol
    
    if c <= 0: return -np.inf
    
    # a. end-of-period assets
    a = m - c
    
    # b. next-period cash-on
    m_plus = par.R*a + par.xi
    
    # c. continuation value
    v_plus = np.zeros(m_plus.size)
    linear_interp.interp_1d_vec(par.grid_m,sol.v[t+1,:],m_plus,v_plus)
            
    # d. value-of-choice
    total = c**(1-par.rho)/(1-par.rho) + par.beta*np.sum(par.xi_w*v_plus)
    return total    
    
def solve(model):
    """ solve model with vfi """
    
    par = model.par
    sol = model.sol
    
    t0 = time.time()
    
    # a. last-period
    sol.c[-1,:] = par.grid_m
    sol.v[-1,:]= par.grid_m**(1-par.rho)/(1-par.rho) 
    
    # b. backwards induction
    for t in reversed(range(par.T-1)):
        for i_m,m in enumerate(par.grid_m):

            obj = lambda c: -value_of_choice(c,t,m,model)
            result = optimize.minimize_scalar(obj,method='bounded',bounds=(0,m))

            sol.c[t,i_m] = result.x
            sol.v[t,i_m]= -result.fun 
    
    print(f'model solved in {elapsed(t0)}')<jupyter_output><empty_output><jupyter_text>The **simulate** function:<jupyter_code>def simulate(model,seed=1917):
    """ simulate model """
    
    par = model.par
    sol = model.sol
    sim = model.sim
    
    t0 = time.time()
    
    # a. shocks
    if not seed is None: np.random.seed(seed)
    sim.xi[:,:] = np.exp(np.random.normal(-0.5*par.sigma_xi**2,par.sigma_xi,size=(par.T,par.simN)))
    
    # a. initial cash-on-hand
    sim.m[0,:] = sim.xi[0,:] 
    
    # c. time loop
    for t in range(par.T):
        
        # i. consumption
        linear_interp.interp_1d_vec(par.grid_m,sol.c[t],sim.m[t,:],sim.c[t,:])

        # ii. end-of-period assets
        sim.a[t,:] = sim.m[t,:]-sim.c[t,:]

        # iii. next-period
        if t < par.T-1: sim.m[t+1] = par.R*sim.a[t,:] + sim.xi[t+1,:]
            
    print(f'model simulated in {elapsed(t0)}')<jupyter_output><empty_output><jupyter_text>Class **Define class:**<jupyter_code>class ConsumptionSavingModelClass(ModelClass):
    
    settings = settings
    setup = setup
    allocate = allocate
    solve = solve
    simulate = simulate<jupyter_output><empty_output><jupyter_text>**Construct instance:**<jupyter_code>model = ConsumptionSavingModelClass(name='myfirstconsavmodel')<jupyter_output><empty_output><jupyter_text>**Check the type inference:**<jupyter_code>print(model)<jupyter_output>Modelclass: ConsumptionSavingModelClass
Name: myfirstconsavmodel

namespaces: ['sol', 'sim', 'par']
other_attrs: []
savefolder: saved
not_floats: ['T', 'Nxi', 'Nm', 'simN']

sol:
 c = ndarray with shape = (100, 100) [dtype: float64]
 v = ndarray with shape = (100, 100) [dtype: float64]
 memory, gb: 0.0

sim:
 m = ndarray with shape = (100, 10000) [dtype: float64]
 c = ndarray with shape = (100, 10000) [dtype: float64]
 a = ndarray with shape = (100, 10000) [dtype: float64]
 xi = ndarray with shape = (100, 10000) [dtype: float64]
 memory, gb: 0.0

par:
 T = 100 [int]
 rho = 2.0 [float]
 beta = 0.96 [float]
 sigma_xi = 0.1 [float]
 Nxi = 8 [int]
 R = 1.04 [float]
 Nm = 100 [int]
 m_min = 1e-06 [float]
 m_max = 20.0 [float]
 m_phi = 1.3 [float]
 simN = 10000 [int]
 grid_m = ndarray with shape = (100,) [dtype: float64]
 xi = ndarray with shape = (8,) [dtype: float64]
 xi_w = ndarray with shape = (8,) [dtype: float64]
 memory, gb: 0.0<jupyter_text>> **Note:** When initializing the model with `ConsumptionSainvgModelClass()` the `.settings()`, `.setup()` and `.allocate()` methods are called, and the types of all variables in `.par`, `.sol` and `.sim` are inferred. These types and the number of dimensions of numpy arrays must not be changed afterwards. **Solve:**<jupyter_code>model.solve()<jupyter_output>model solved in 4.6 secs<jupyter_text>**Simulate:**<jupyter_code>model.simulate()
par = model.par
sol = model.sol
sim = model.sim<jupyter_output><empty_output><jupyter_text>**Consumption functions:**<jupyter_code>fig = plt.figure()
ax = fig.add_subplot(1,1,1)

for t in [par.T-1, par.T-2, par.T-6, par.T-11, 50, 25, 0]:
    if t >= 0: ax.plot(par.grid_m,sol.c[t,:],'-o',ms=4,label=f'$t = {t}$')

# layout
ax.legend(frameon=True)
ax.set_xlim([0,5])
ax.set_ylim([0,5])
ax.set_xlabel('$m_t$')
ax.set_ylabel('$c_t$');<jupyter_output><empty_output><jupyter_text>**Time profile of end-of-period assets:**<jupyter_code>par = model.par
sol = model.sol

fig = plt.figure()
ax = fig.add_subplot(1,1,1)

ax.plot(np.arange(par.T),np.mean(sim.a,axis=1))

# layout
ax.set_xlabel('time')
ax.set_ylabel('$a_t$');<jupyter_output><empty_output><jupyter_text>**Histograms of end-of-period assets:**<jupyter_code>fig = plt.figure()
ax = fig.add_subplot(1,1,1)

for t in [25,50,75]:
    if t < par.T-1: ax.hist(sim.a[t,:],bins=100,alpha=0.5,label=f'$t = {t}$',density=True)

# layout
ax.legend(frameon=True)
ax.set_xlabel('$m_t$');<jupyter_output><empty_output><jupyter_text>Inherited methods **1. update at initilization:**<jupyter_code>low_beta = ConsumptionSavingModelClass(name='lowbeta',par={'beta':0.9})
print(low_beta.par.beta)<jupyter_output>0.9<jupyter_text>**Note:** Happens after `.setup` is called, but before `.allocate` is called **2. save/load:**<jupyter_code>model.save() # see saved\myfirstconsavmodel*
print(f'checksum: {np.sum(model.sol.c):.8f}')
model_loaded = ConsumptionSavingModelClass(name='myfirstconsavmodel',load=True)
print(f'checksum: {np.sum(model_loaded.sol.c):.8f}')<jupyter_output>checksum: 12104.94231004<jupyter_text>**3. copy**<jupyter_code>model_copy = model_loaded.copy()
print(f'name: {model_copy.name}')
print(f'checksum: {np.sum(model_copy.sol.c):.8f}')<jupyter_output>name: myfirstconsavmodel_copy
checksum: 12104.94231004<jupyter_text>**4. as dictionary** (can be pickled)<jupyter_code>model_dict = model.as_dict()
model_fromdict = ConsumptionSavingModelClass(name='fromdict',from_dict=model_dict)
print(f'name: {model_fromdict.name}')
print(f'checksum: {np.sum(model_fromdict.sol.c):.8f}')<jupyter_output>name: fromdict
checksum: 12104.94231004<jupyter_text>Numba The beauty of the **ModelClass** from **ConSav** is that it is [Numba]((http://numba.pydata.org/)) prepared, so it is very easy to call numba jitted functions. Numba automatically convert Python and NumPy code into fast machine code. Specifically, this is done for all functions preceded (decorated) with ``@njit`` or ``@njit(parallel=True)`` when parallization is used.**Note:** For more info about **Numba** see [here](https://numeconcopenhagen.netlify.com/lectures/The_need_for_speedNumba). **New functions:**<jupyter_code>@njit
def value_of_choice_jit(c,t,m,sol,par):
            
    # a. end-of-period assets
    a = m - c
    
    # b. next-period cash-on
    m_plus = par.R*a + par.xi
    
    # c. continuation value
    v_plus = np.zeros(m_plus.size)
    linear_interp.interp_1d_vec(par.grid_m,sol.v[t+1,:],m_plus,v_plus)
            
    # d. value-of-choice
    total = c**(1-par.rho)/(1-par.rho) + par.beta*np.sum(par.xi_w*v_plus)
    return -total    

@njit(parallel=True)
def solve_(t,sol,par):
    
    c = sol.c[t]
    v = sol.v[t]
    
    for i_m in prange(par.Nm):
        m = par.grid_m[i_m]
        c[i_m] = golden_section_search.optimizer(value_of_choice_jit,1e-8,m,args=(t,m,sol,par))
        v[i_m] = -value_of_choice_jit(c[i_m],t,m,sol,par)
        
def solve_jit(model):
    
    par = model.par
    sol = model.sol
    
    t0 = time.time()
    
    # a. last-period
    sol.c[-1,:] = par.grid_m
    sol.v[-1,:]= par.grid_m**(1-par.rho)/(1-par.rho) 
    
    # b. backwards induction
    with jit(model) as model_jit: 
        for t in reversed(range(par.T-1)):
            solve_(t,model_jit.sol,model_jit.par)
        
    print(f'model solved in {elapsed(t0)}')<jupyter_output><empty_output><jupyter_text>**New model class:**<jupyter_code>class ConsumptionSavingModelClass_jit(ModelClass):
    
    settings = settings
    setup = setup
    allocate = allocate
    solve = solve_jit
    simulate = simulate<jupyter_output><empty_output><jupyter_text>**Solve the model with numba:**1. On the first run a few seconds is spend on analyzing the code2. From the second run onwards it is really fast<jupyter_code>t0 = time.time()
model_jit = ConsumptionSavingModelClass_jit(name='jit')
print(f'model setup done in {elapsed(t0)}')

for _ in range(3):
    model_jit.solve()<jupyter_output>model setup done in 0.0 secs
model solved in 2.6 secs
model solved in 0.2 secs
model solved in 0.2 secs<jupyter_text>**Note:** When copying a model, the code is NOT re-analyzed<jupyter_code>model_jit_alt = model_jit.copy()

for _ in range(3):
    model_jit_alt.solve()<jupyter_output>model solved in 0.2 secs
model solved in 0.2 secs
model solved in 0.2 secs<jupyter_text>**Compare solutions:**<jupyter_code>fig = plt.figure()
ax = fig.add_subplot(1,1,1)

ax.plot(par.grid_m,model.sol.c[0,:],'-o',ms=4,label=f'pure python')
ax.plot(par.grid_m,model_jit.sol.c[0,:],'-o',ms=4,label=f'with numba')

# layout
ax.legend(frameon=True)
ax.set_xlim([0,5])
ax.set_ylim([0,5])
ax.set_xlabel('$m_t$')
ax.set_ylabel('$c_t$');<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15850080680651313, 'fraction_numerical': 0.023764119113979756, 'mean_word_length': 3.987198244330651, 'pattern_counts': {'":': 0, '<': 88, '<?xml version=': 0, '>': 87, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '23', 'dir': 'jupyter-structured-clean-dedup', 'id': '355690', 'license': 'MIT', 'n_tokens_mistral': 5497, 'n_tokens_neox': 5231, 'n_words': 1114, 'path': '00. DynamicProgramming/02. Your First Consumption-Saving Model.ipynb', 'repo_name': 'alanlujan91/ConsumptionSavingNotebooks', 'size': '13634'}