Model solution#

In this exercise, we will go through the model solution.

Warning

Remember that ams.func is the sintax that calls the file with the functions see before.

Packages#

First, calls numpy and my functions file so we can use it in our model solution.

import numpy as np
import ams_functions as ams

Parameters#

To start our solution, we need to determine the parameters. The names here can be a little bit confusion, but the age_max indicate how many adding yeas the child can study. This means that we are simulated children that are 6 years old and can study until they are 16.

age_start = 6
age_max   = 10

Creating matrices#

Before the model solution loop, we need to create empty matrices to store our results. I commented the code so you can see what is the matrix about.

# Probabilities of working
probs_w = np.empty((age_max,age_max)) 
probs_w[:] = np.nan         

# Probabilities of school
probs_s = np.empty((age_max,age_max)) 
probs_s[:] = np.nan   

# Threshold matrix
eps_t = np.empty((age_max,age_max)) 
eps_t[:] = np.nan         

# Expected values matrix for L = 1 
eve_w  = np.empty((age_max,age_max)) 
eve_w[:] = np.nan              

# Expected values matrix for L = 0 
eve_s = np.empty((age_max,age_max)) 
eve_s[:] = np.nan   

# Value function matrix
vf = np.empty((age_max,age_max)) 
vf[:] = np.nan  

Model solution loop#

The model solution is solved backward and follows these steps:

  1. Solve for the threshold that makes the value function for going to school bigger.

  2. Use the result from 1 to obtain the probability of going to school.

  3. Calculate the truncated distribution value for going to school

  4. Calculate the EV (expected values in the period ahead)

  5. Repeat until the end of the loop

In this case, for the last period we will use the terminal value defined as:

(1)#\[\begin{align} V(ed_{i,18}) = \frac{\alpha_{1}}{1 + exp(-\alpha_{2}*ed_{i,18})} \end{align}\]

This gives the value of having an specific level of education in the last period

for age in range(age_max,-1,-1):
    for edu in range(0,age_max):
        if edu >= age:
            pass
        else:       
            # For the threshold we will use the terminal value to capture the value of an specifiy education level in terms of utility  
               
            if age == age_max: #last period situation
                tv_s =  ams.prob_progress(edu)*(ams.terminal_v(edu = edu+1)) + (1 - ams.prob_progress(edu))*(ams.terminal_v(edu = edu)) 
                tv_w =  ams.terminal_v(edu = edu)
                eps_t[age-1,edu] = ams.value_function(age,edu,1, EV = tv_s) - ams.value_function(age,edu,0, EV = tv_w)
                                   
            else:                #Other periods will take the difference between value functions
                tv_s = ams.prob_progress(edu)*(vf[age, edu+1]) + (1 - ams.prob_progress(edu))*(vf[age,edu])
                tv_w = vf[age, edu]
                eps_t[age-1,edu] = ams.value_function(age,edu,1, EV = tv_s) - ams.value_function(age,edu,0, EV = tv_w)
            
            # Now we will calcute the probabilities of schooling and working based on the thresholds                     
            probs_s[age-1,edu] = ams.logistic(eps = eps_t[age-1,edu])
            probs_w[age-1,edu] = 1 - probs_s[age-1,edu]
            
            # Using the thresholds we will calculate the EVs
            if age == age_max:
                u_t = ams.trunc_change(eps = eps_t[age-1,edu])
                eve_s[age-1,edu] = ams.value_function(age,edu,1, EV = tv_s) + ams.trunc_school(ϵ_threshold = eps_t[age-1,edu], u_threshold = u_t)
                eve_w[age-1,edu] = ams.value_function(age,edu,0, EV = tv_w)
            else:
                u_t = ams.trunc_change(eps = eps_t[age-1,edu])
                eve_s[age-1,edu] = ams.value_function(age,edu,1, EV = tv_s) + ams.trunc_school(ϵ_threshold = eps_t[age-1,edu], u_threshold = u_t)
                eve_w[age-1,edu] = ams.value_function(age,edu,0, EV = tv_w)
            
            vf[age-1,edu] = probs_s[age-1,edu]*eve_s[age-1,edu] + probs_w[age-1,edu]*eve_w[age-1,edu]

Let’s some results. The probability of going to school matrix looks like this:

probs_s
array([[0.99999993,        nan,        nan,        nan,        nan,
               nan,        nan,        nan,        nan,        nan],
       [0.99999884, 0.99999981,        nan,        nan,        nan,
               nan,        nan,        nan,        nan,        nan],
       [0.99997971, 0.99999692, 0.99999955,        nan,        nan,
               nan,        nan,        nan,        nan,        nan],
       [0.99962732, 0.99994308, 0.99999232, 0.99999898,        nan,
               nan,        nan,        nan,        nan,        nan],
       [0.99586383, 0.99880728, 0.99984954, 0.99998211, 0.99999649,
               nan,        nan,        nan,        nan,        nan],
       [0.91581254, 0.98644322, 0.9962496 , 0.99962788, 0.9999358 ,
        0.99998729,        nan,        nan,        nan,        nan],
       [0.65754572, 0.76395663, 0.96031838, 0.98826779, 0.99857002,
        0.99975688, 0.99994754,        nan,        nan,        nan],
       [0.33117145, 0.39518156, 0.55636764, 0.91813627, 0.94320286,
        0.99389726, 0.9989383 , 0.99971736,        nan,        nan],
       [0.13441563, 0.15878852, 0.11884147, 0.7362356 , 0.67868231,
        0.76872596, 0.97157708, 0.99372576, 0.99781425,        nan],
       [0.05324222, 0.05684499, 0.03267233, 0.43056996, 0.28367671,
        0.27255018, 0.50621805, 0.88952229, 0.94513499, 0.97851892]])

The EV matrix is:

vf
array([[154.2996938 ,          nan,          nan,          nan,
                 nan,          nan,          nan,          nan,
                 nan,          nan],
       [147.40631343, 165.27423062,          nan,          nan,
                 nan,          nan,          nan,          nan,
                 nan,          nan],
       [142.41339016, 160.14706687, 179.02221917,          nan,
                 nan,          nan,          nan,          nan,
                 nan,          nan],
       [139.62337579, 157.12161177, 175.85698874, 195.86431707,
                 nan,          nan,          nan,          nan,
                 nan,          nan],
       [139.31015202, 156.54370737, 175.01510424, 194.87866649,
        211.0367933 ,          nan,          nan,          nan,
                 nan,          nan],
       [141.18076829, 158.80179566, 176.88834836, 196.46108271,
        212.47273586, 228.6402713 ,          nan,          nan,
                 nan,          nan],
       [146.3594149 , 163.55689397, 182.03683495, 201.05416816,
        216.74777283, 232.76714143, 248.15710453,          nan,
                 nan,          nan],
       [153.3799503 , 171.87520601, 189.96149747, 209.48071548,
        224.37003875, 240.03721847, 255.27493485, 268.60861758,
                 nan,          nan],
       [161.33329164, 181.63603512, 201.23849672, 220.60613972,
        236.52013511, 251.09686622, 265.88547942, 279.05383333,
        289.06793576,          nan],
       [169.66802391, 192.06435084, 213.670688  , 233.8340728 ,
        251.59177304, 267.19573071, 280.76515423, 293.43892839,
        303.18467998, 311.39752683]])

Using our function file, you can use the following code to run the loop we just saw

probs_w, probs_s, eps_t, eve_w, eve_s, vf = ams.solution()