From 70b271712695f80caa0fb69dfe9f146e969084a1 Mon Sep 17 00:00:00 2001 From: Richard Evans Date: Wed, 23 Oct 2024 00:20:06 +0800 Subject: [PATCH 1/4] Added solutions to chapter problems --- .../content/UNtutorial/getting_started.md | 2 +- .../UNtutorial/solutions/3perOG/README.md | 14 + .../content/UNtutorial/solutions/3perOG/SS.py | 628 +++++++++++++++++ .../UNtutorial/solutions/3perOG/TPI.py | 640 ++++++++++++++++++ .../UNtutorial/solutions/3perOG/execute.py | 200 ++++++ 5 files changed, 1483 insertions(+), 1 deletion(-) create mode 100644 docs/book/content/UNtutorial/solutions/3perOG/README.md create mode 100644 docs/book/content/UNtutorial/solutions/3perOG/SS.py create mode 100644 docs/book/content/UNtutorial/solutions/3perOG/TPI.py create mode 100755 docs/book/content/UNtutorial/solutions/3perOG/execute.py diff --git a/docs/book/content/UNtutorial/getting_started.md b/docs/book/content/UNtutorial/getting_started.md index c70b1cd..33e5baa 100644 --- a/docs/book/content/UNtutorial/getting_started.md +++ b/docs/book/content/UNtutorial/getting_started.md @@ -32,7 +32,7 @@ For the October 21-25, 2024 United Nations `OG-PHL` training in Manila, we will - Theory: ["Simple" 3-period-lived agent model](https://eapd-drb.github.io/OG-PHL/content/UNtutorial/3perOG.html) * - Tue. - Morning - - Review 3-period-lived-agent exercises
Review OG-Core and OG-PHL modules
Quick Git and GitHub workflow review + - Review 3-period-lived-agent exercises (solutions in [this folder](https://github.com/EAPD-DRB/OG-PHL/blob/main/docs/book/content/UNtutorial/Solutions/3perOG/))
Review OG-Core and OG-PHL modules
Quick Git and GitHub workflow review * - - Afternoon - Running OG-PHL, inputs, outputs
Calibrating OG-PHL, current state, still to do diff --git a/docs/book/content/UNtutorial/solutions/3perOG/README.md b/docs/book/content/UNtutorial/solutions/3perOG/README.md new file mode 100644 index 0000000..def5200 --- /dev/null +++ b/docs/book/content/UNtutorial/solutions/3perOG/README.md @@ -0,0 +1,14 @@ +# Solutions to the exercises in the 3-period-lived OG model chapter in the UN Tutorial section +This folder contains the python scripts that compute the solutions to [Exercise 1](https://eapd-drb.github.io/OG-PHL/content/UNtutorial/3perOG.html#ExerUNtut_3perOGfeas), [Exercise 2](https://eapd-drb.github.io/OG-PHL/content/UNtutorial/3perOG.html#ExerUNtut_3perOGSS), and [Exercise 3](https://eapd-drb.github.io/OG-PHL/content/UNtutorial/3perOG.html#ExerUNtut_3perOGTPI) in the "[A 'Simple' 3-period-lived agent OG model](https://eapd-drb.github.io/OG-PHL/content/UNtutorial/3perOG.html)". In addition to this `README.md` file, this folder contains the following three Python scripts. +- [`execute.py`](https://github.com/EAPD-DRB/OG-PHL/blob/main/docs/book/content/UNtutorial/solutions/3perOG/execute.py). This script runs everything and delivers the solutions to the exercises. It calls the following two modules [`SS.py`](https://github.com/EAPD-DRB/OG-PHL/blob/main/docs/book/content/UNtutorial/solutions/3perOG/SS.py) and [`TPI.py`](https://github.com/EAPD-DRB/OG-PHL/blob/main/docs/book/content/UNtutorial/solutions/3perOG/TPI.py). +- [`SS.py`](https://github.com/EAPD-DRB/OG-PHL/blob/main/docs/book/content/UNtutorial/solutions/3perOG/SS.py). Python module that has the functions for solving the feasibility questions of [Exercise 1](https://eapd-drb.github.io/OG-PHL/content/UNtutorial/3perOG.html#ExerUNtut_3perOGfeas) and the steady-state equilibrium questions of [Exercise 2](https://eapd-drb.github.io/OG-PHL/content/UNtutorial/3perOG.html#ExerUNtut_3perOGSS). +- [`TPI.py`](https://github.com/EAPD-DRB/OG-PHL/blob/main/docs/book/content/UNtutorial/solutions/3perOG/TPI.py). Python module that has the functions for solving the transition path equilibrium questions of [Exercise 3](https://eapd-drb.github.io/OG-PHL/content/UNtutorial/3perOG.html#ExerUNtut_3perOGTPI). + +You can run these computations that solve for the solutions to the three exercises by doing the following steps. +1. Open your terminal windown (Mac or Linux) or Anaconda prompt (Windows). +2. Activate your `ogphl-dev` conda environment: `conda activate ogphl-dev` +3. Navigate to your `OG-PHL` repository directory. +4. Navigate to the `/docs/book/content/UNtutorial/solutions/3perOG/` subfolder of this repository. +5. Run the `execute.py` script: `python execute.py`. + +Running the `execute.py` script will print all the answers to the exercises in the output of the terminal. It will also create an `images` folder in this directory that has some of the images from the steady state and transition path equilibrium solutions. diff --git a/docs/book/content/UNtutorial/solutions/3perOG/SS.py b/docs/book/content/UNtutorial/solutions/3perOG/SS.py new file mode 100644 index 0000000..231bd87 --- /dev/null +++ b/docs/book/content/UNtutorial/solutions/3perOG/SS.py @@ -0,0 +1,628 @@ +''' +------------------------------------------------------------------------ +This module contains the functions used to solve the steady state for +the model with 3-period lived agents and exogenous labor from Chapter 5 +of the OG textbook. + +This Python module calls the following function(s): + print_time() + get_cvec() + get_L() + get_K() + get_w() + get_r() + get_Y() + get_C() + feasible() + EulerSys() + get_b_errors() + get_SS() +------------------------------------------------------------------------ +''' +# Import packages +import time +import numpy as np +import scipy.optimize as opt +import matplotlib.pyplot as plt +from matplotlib.ticker import MultipleLocator +import os + + +# Define functions + + +def print_time(seconds, type): + ''' + -------------------------------------------------------------------- + Takes a total amount of time in seconds and prints it in terms of + more readable units (days, hours, minutes, seconds) + -------------------------------------------------------------------- + INPUTS: + seconds = scalar > 0, total amount of seconds + type = string, either "SS" or "TPI" + + OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION: + + OBJECTS CREATED WITHIN FUNCTION: + secs = scalar > 0, remainder number of seconds + mins = integer >= 1, remainder number of minutes + hrs = integer >= 1, remainder number of hours + days = integer >= 1, number of days + + FILES CREATED BY THIS FUNCTION: None + + RETURNS: Nothing + -------------------------------------------------------------------- + ''' + if seconds < 60: # seconds + secs = round(seconds, 4) + print(type + ' computation time: ' + str(secs) + ' sec') + elif seconds >= 60 and seconds < 3600: # minutes + mins = int(seconds / 60) + secs = round(((seconds / 60) - mins) * 60, 1) + print(type + ' computation time: ' + str(mins) + ' min, ' + + str(secs) + ' sec') + elif seconds >= 3600 and seconds < 86400: # hours + hrs = int(seconds / 3600) + mins = int(((seconds / 3600) - hrs) * 60) + secs = round(((seconds / 60) - hrs * 60 - mins) * 60, 1) + print(type + ' computation time: ' + str(hrs) + ' hrs, ' + + str(mins) + ' min, ' + str(secs) + ' sec') + elif seconds >= 86400: # days + days = int(seconds / 86400) + hrs = int(((seconds / 86400) - days) * 24) + mins = int(((seconds / 3600) - days * 24 - hrs) * 60) + secs = round( + ((seconds / 60) - days * 24 * 60 - hrs * 60 - mins) * 60, 1) + print(type + ' computation time: ' + str(days) + ' days, ' + + str(hrs) + ' hrs, ' + str(mins) + ' min, ' + + str(secs) + ' sec') + + +def get_cvec(r, w, bvec, nvec): + ''' + -------------------------------------------------------------------- + Generates vector of lifetime steady-state consumptions given savings + decisions, parameters, and the corresponding steady-state interest + rate and wage. + -------------------------------------------------------------------- + INPUTS: + r = scalar > 0, steady-state interest rate + w = scalar > 0, steady-state wage + bvec = (S-1,) vector, distribution of savings b_{s+1} + nvec = (S,) vector, exogenous labor supply n_{s} + + OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION: + + OBJECTS CREATED WITHIN FUNCTION: + b_s = (S,) vector, b_init in first element and bvec in last S-1 + elements + b_sp1 = (S,) vector, bvec in first S-1 elements and 0 in last + element + cvec = (S,) vector, consumption by age c_s + c_cnstr = (S,) Boolean vector, =True if c_s <= 0 + + FILES CREATED BY THIS FUNCTION: None + + RETURNS: cvec, c_cnstr + -------------------------------------------------------------------- + ''' + b_s = np.append([0], bvec) + b_sp1 = np.append(bvec, [0]) + cvec = (1 + r) * b_s + w * nvec - b_sp1 + if cvec.min() <= 0: + print('get_cvec() warning: distribution of savings and/or ' + + 'parameters created c<=0 for some agent(s)') + c_cnstr = cvec <= 0 + + return cvec, c_cnstr + + +def get_L(nvec): + ''' + -------------------------------------------------------------------- + Solve for aggregate labor L + -------------------------------------------------------------------- + INPUTS: + nvec = (3,) vector, exogenous labor supply values (n_1, n_2, n_3) + + OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION: None + + OBJECTS CREATED WITHIN FUNCTION: + L = scalar > 0, aggregate labor + + FILES CREATED BY THIS FUNCTION: None + + RETURNS: L + -------------------------------------------------------------------- + ''' + L = nvec.sum() + + return L + + +def get_K(barr): + ''' + -------------------------------------------------------------------- + Solve for steady-state aggregate capital stock K or time path of + aggregate capital stock K_t + -------------------------------------------------------------------- + INPUTS: + barr = (2,) vector or (2, T+S-2) matrix, values for steady-state + savings (b_2, b_3) or time path of distribution of savings + + OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION: None + + OBJECTS CREATED WITHIN FUNCTION: + K = scalar or (T+S-2,) vector, steady-state aggregate capital + stock or time path of aggregate capital stock + K_cnstr = Boolean or (T+S-2) Boolean, =True if K <= 0 or if K_t <= 0 + + FILES CREATED BY THIS FUNCTION: None + + RETURNS: K, K_cnstr + -------------------------------------------------------------------- + ''' + if barr.ndim == 1: # This is the steady-state case + K = barr.sum() + K_cnstr = K <= 0 + if K_cnstr: + print('get_K() warning: distribution of savings and/or ' + + 'parameters created K<=0 for some agent(s)') + + elif barr.ndim == 2: # This is the time path case + K = barr.sum(axis=0) + K_cnstr = K <= 0 + if K.min() <= 0: + print('Aggregate capital constraint is violated K<=0 for ' + + 'some period in time path.') + + return K, K_cnstr + + +def get_w(params, K, L): + ''' + -------------------------------------------------------------------- + Solve for steady-state wage w or time path of wages w_t + -------------------------------------------------------------------- + INPUTS: + params = length 2 tuple, (A, alpha) + A = scalar > 0, total factor productivity + alpha = scalar in (0, 1), capital share of income + K = scalar > 0 or (T+S-2,) vector, steady-state aggregate + capital stock or time path of the aggregate capital stock + L = scalar > 0 or (T+S-2,) vector, steady-state aggregate + labor or time path of aggregate labor + + OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION: None + + OBJECTS CREATED WITHIN FUNCTION: + w = scalar > 0 or (T+S-2) vector, steady-state wage or time path of + wage + + FILES CREATED BY THIS FUNCTION: None + + RETURNS: w + -------------------------------------------------------------------- + ''' + A, alpha = params + w = (1 - alpha) * A * ((K / L) ** alpha) + + return w + + +def get_r(params, K, L): + ''' + -------------------------------------------------------------------- + Solve for steady-state interest rate r or time path of interest + rates r_t + -------------------------------------------------------------------- + INPUTS: + params = length 3 tuple, (A, alpha, delta) + A = scalar > 0, total factor productivity + alpha = scalar in (0, 1), capital share of income + delta = scalar in (0, 1), per period depreciation rate + K = scalar > 0 or (T+S-2,) vector, steady-state aggregate + capital stock or time path of the aggregate capital stock + L = scalar > 0 or (T+S-2,) vector, steady-state aggregate + labor or time path of aggregate labor + + OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION: None + + OBJECTS CREATED WITHIN FUNCTION: + r = scalar > 0 or (T+S-2) vector, steady-state interest rate or time + path of interest rate + + FILES CREATED BY THIS FUNCTION: None + + RETURNS: r + -------------------------------------------------------------------- + ''' + A, alpha, delta = params + r = alpha * A * ((L / K) ** (1 - alpha)) - delta + + return r + + +def get_Y(params, K, L): + ''' + -------------------------------------------------------------------- + Solve for steady-state aggregate output Y or time path of aggregate + output Y_t + -------------------------------------------------------------------- + INPUTS: + params = length 2 tuple, production function parameters + (A, alpha) + A = scalar > 0, total factor productivity + alpha = scalar in (0,1), capital share of income + K = scalar > 0 or (T+S-2,) vector, aggregate capital stock + or time path of the aggregate capital stock + L = scalar > 0 or (T+S-2,) vector, aggregate labor or time + path of the aggregate labor + + OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION: None + + OBJECTS CREATED WITHIN FUNCTION: + Y = scalar > 0 or (T+S-2,) vector, aggregate output (GDP) or + time path of aggregate output (GDP) + + FILES CREATED BY THIS FUNCTION: None + + RETURNS: Y + -------------------------------------------------------------------- + ''' + A, alpha = params + Y = A * (K ** alpha) * (L ** (1 - alpha)) + + return Y + + +def get_C(carr): + ''' + -------------------------------------------------------------------- + Solve for steady-state aggregate consumption C or time path of + aggregate consumption C_t + -------------------------------------------------------------------- + INPUTS: + carr = (S,) vector or (S, T) matrix, distribution of consumption c_s + in steady state or time path for the distribution of + consumption + + OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION: None + + OBJECTS CREATED WITHIN FUNCTION: + C = scalar > 0 or (T,) vector, aggregate consumption or time path of + aggregate consumption + + Returns: C + -------------------------------------------------------------------- + ''' + if carr.ndim == 1: + C = carr.sum() + elif carr.ndim == 2: + C = carr.sum(axis=0) + + return C + + +def feasible(params, bvec): + ''' + -------------------------------------------------------------------- + Check whether a vector of steady-state savings is feasible in that + it satisfies the nonnegativity constraints on consumption in every + period c_s > 0 and that the aggregate capital stock is strictly + positive K > 0 + -------------------------------------------------------------------- + INPUTS: + params = length 4 tuple, (nvec, A, alpha, delta) + nvec = (3,) vector, exogenous labor supply values (n_1, n_2, n_3) + A = scalar > 0, total factor productivity + alpha = scalar in (0, 1), capital share of income + delta = scalar in (0, 1), per period depreciation rate + bvec = (2,) vector, values of steady-state savings ([b_2, b_3]) + + OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION: + get_L() + get_K() + get_w() + get_r() + get_cvec() + + OBJECTS CREATED WITHIN FUNCTION: + L = scalar > 0, steady-state aggregate labor + K = scalar, steady-state aggregate capital stock + K_cnstr = Boolean, =True if K <= 0 + w_params = length 2 tuple, (A, alpha) + w = scalar, steady-state wage + r_params = length 3 tuple, (A, alpha, delta) + r = scalar, steady-state interest rate + cvec = (3,) vector, steady-state consumption by age + c_cnstr = (3,) Boolean vector, =True for elements for which c_s<=0 + b_cnstr = (2,) Boolean, =True for elements for which b_s causes a + violation of the nonnegative consumption constraint. + + FILES CREATED BY THIS FUNCTION: None + + RETURNS: b_cnstr, c_cnstr, K_cnstr + -------------------------------------------------------------------- + ''' + nvec, A, alpha, delta = params + L = get_L(nvec) + K, K_cnstr = get_K(bvec) + if not K_cnstr: + w_params = (A, alpha) + w = get_w(w_params, K, L) + r_params = (A, alpha, delta) + r = get_r(r_params, K, L) + cvec, c_cnstr = get_cvec(r, w, bvec, nvec) + b_cnstr = c_cnstr[:-1] + c_cnstr[1:] + else: + c_cnstr = np.ones(cvec.shape[0], dtype=bool) + b_cnstr = np.ones(cvec.shape[0] - 1, dtype=bool) + + return b_cnstr, c_cnstr, K_cnstr + + +def EulerSys(bvec, *args): + ''' + -------------------------------------------------------------------- + Generates vector of all Euler errors that characterize optimal + lifetime decisions + -------------------------------------------------------------------- + INPUTS: + bvec = (2,) vector, distribution of savings b_{s+1} + args = length 8 tuple, + (beta, sigma, nvec, L, A, alpha, delta, EulDiff) + beta = scalar in [0,1), discount factor + sigma = scalar > 0, coefficient of relative risk aversion + nvec = (S,) vector, exogenous labor supply n_{s} + L = scalar > 0, exogenous aggregate labor + A = scalar > 0, total factor productivity parameter in firms' + production function + alpha = scalar in (0,1), capital share of income + delta = scalar in [0,1], model-period depreciation rate of capital + EulDiff = Boolean, =True if want difference version of Euler errors + beta*(1+r)*u'(c2) - u'(c1), =False if want ratio version + [beta*(1+r)*u'(c2)]/[u'(c1)] - 1 + + OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION: + get_K() + get_r() + get_w() + get_cvec() + get_b_errors() + + OBJECTS CREATED WITHIN FUNCTION: + K = scalar > 0, aggregate capital stock + K_constr = Boolean, =True if K<=0 for given bvec + b_err_vec = (S-1,) vector, vector of Euler errors + r_params = length 3 tuple, (A, alpha, delta) + r = scalar > 0, interest rate + w_params = length 2 tuple, (A, alpha) + w = scalar > 0, wage + cvec = (S,) vector, consumption c_s for each age-s agent + c_constr = (S,) Boolean vector, =True if c<=0 for given bvec + b_err_params = length 2 tuple, (beta, sigma) + + FILES CREATED BY THIS FUNCTION: None + + RETURNS: b_errors + -------------------------------------------------------------------- + ''' + beta, sigma, nvec, L, A, alpha, delta, EulDiff = args + K, K_cnstr = get_K(bvec) + if K_cnstr: + b_err_vec = 1000. * np.ones(nvec.shape[0] - 1) + else: + r_params = (A, alpha, delta) + r = get_r(r_params, K, L) + w_params = (A, alpha) + w = get_w(w_params, K, L) + cvec, c_cnstr = get_cvec(r, w, bvec, nvec) + b_err_params = (beta, sigma) + b_err_vec = get_b_errors(b_err_params, r, cvec, c_cnstr, + EulDiff) + + return b_err_vec + + +def get_b_errors(params, r, cvec, c_cnstr, diff): + ''' + -------------------------------------------------------------------- + Generates vector of dynamic Euler errors that characterize the + optimal lifetime savings decision. Because this function is used for + solving for lifetime decisions in both the steady-state and in the + transition path, lifetimes will be of varying length. Lifetimes in + the steady-state will be S periods. Lifetimes in the transition path + will be p in [2, S] periods + -------------------------------------------------------------------- + INPUTS: + params = length 2 tuple, (beta, sigma) + beta = scalar in (0,1), discount factor + sigma = scalar > 0, coefficient of relative risk aversion + r = scalar > 0 or (p,) vector, steady-state interest rate or + time path of interest rates + cvec = (p,) vector, distribution of consumption by age c_p + c_cnstr = (p,) Boolean vector, =True if c<=0 for given bvec + diff = boolean, =True if use simple difference Euler + errors. Use percent difference errors otherwise. + + OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION: None + + OBJECTS CREATED WITHIN FUNCTION: + mu_c = (p-1,) vector, marginal utility of current consumption + mu_cp1 = (p-1,) vector, marginal utility of next period consumpt'n + b_errors = (p-1,) vector, Euler errors characterizing optimal + savings bvec + + FILES CREATED BY THIS FUNCTION: None + + RETURNS: b_errors + -------------------------------------------------------------------- + ''' + beta, sigma = params + # Make each negative consumption artifically positive + cvec[c_cnstr] = 9999. + mu_c = cvec[:-1] ** (-sigma) + mu_cp1 = cvec[1:] ** (-sigma) + if diff: + b_errors = (beta * (1 + r) * mu_cp1) - mu_c + b_errors[c_cnstr[:-1]] = 9999. + b_errors[c_cnstr[1:]] = 9999. + else: + b_errors = ((beta * (1 + r) * mu_cp1) / mu_c) - 1 + b_errors[c_cnstr[:-1]] = 9999. / 100 + b_errors[c_cnstr[1:]] = 9999. / 100 + + return b_errors + + +def get_SS(params, bvec_guess, graphs): + ''' + -------------------------------------------------------------------- + Solve for the steady-state solution of the 3-period-lived agent OG + model with exogenous labor supply + -------------------------------------------------------------------- + INPUTS: + params = length 9 tuple, (beta, sigma, nvec, L, A, alpha, delta, + SS_tol, EulDiff) + beta = scalar in (0,1), discount factor for each model period + sigma = scalar > 0, coefficient of relative risk aversion + nvec = [S,] vector, exogenous labor supply n_{s,t} + L = scalar > 0, exogenous aggregate labor + A = scalar > 0, total factor productivity parameter in + firms' production function + alpha = scalar in (0,1), capital share of income + delta = scalar in [0,1], model-period depreciation rate of + capital + SS_tol = scalar > 0, tolerance level for steady-state fsolve + EulDiff = Boolean, =True if want difference version of Euler + errors beta*(1+r)*u'(c2) - u'(c1), =False if want ratio + version [beta*(1+r)*u'(c2)]/[u'(c1)] - 1 + bvec_guess = (2,) vector, initial guesses for steady-state savings + ([b_2, b_3]) + graphs = Boolean, =True if output steady-state graphs + + OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION: + + OBJECTS CREATED WITHIN FUNCTION: + f_params = length 4 tuple, (nvec, A, alpha, delta) + b1_cnstr = (2,) Boolean vector, =True if b_s causes negative + consumption c_s <= 0 or negative aggregate capital + stock K <= 0 + c1_cnstr = (3,) Boolean vector, =True for elements of negative + consumption c_s <= 0 + K1_cnstr = Boolean, =True if K <= 0 + eul_args = length 8 tuple, (beta, sigma, nvec, L, A, alpha, + delta, EulDiff) + b_ss = (2,) vector, steady-state distribution of savings + K_ss = scalar > 0, steady-state aggregate capital stock + K_cnstr = Boolean, =True if K_ss <= 0 + r_params = length 3 tuple, (A, alpha, delta) + r_ss = scalar > 0, steady-state interest rate + w_params = length 2 tuple, (A, alpha) + w_ss = scalar > 0, steady-state wage + c_ss = (3,) vector, steady-state distribution of consumption + c_cnstr = (3,) Boolean vector, + Y_params = length 2 tuple, (A, alpha) + Y_ss = scalar > 0, steady-state aggregate output (GDP) + C_ss = scalar > 0, steady-state aggregate consumption + b_err_params = length 2 tuple, (beta, sigma) + EulErr_ss = (2,) vector, steady-state Euler errors + RCerr_ss = scalar, resource constraint error Y_t - C_t - I_t + ss_time = scalar > 0, time elapsed in seconds + ss_output = length 10 dictionary, {b_ss, c_ss, w_ss, r_ss, K_ss, + Y_ss, C_ss, EulErr_ss, RCerr_ss, ss_time} + + FILES CREATED BY THIS FUNCTION: + SS_bc.png + + RETURNS: ss_output + -------------------------------------------------------------------- + ''' + start_time = time.time() + beta, sigma, nvec, L, A, alpha, delta, SS_tol, EulDiff = params + f_params = (nvec, A, alpha, delta) + b1_cnstr, c1_cnstr, K1_cnstr = feasible(f_params, bvec_guess) + if K1_cnstr is True or c1_cnstr.max() is True: + err = ("Initial guess problem: " + + "One or more constraints not satisfied.") + print("K1_cnstr: ", K1_cnstr) + print("c1_cnstr: ", c1_cnstr) + raise RuntimeError(err) + else: + eul_args = (beta, sigma, nvec, L, A, alpha, delta, EulDiff) + results_bss = opt.root(EulerSys, bvec_guess, args=(eul_args), + tol=SS_tol) + b_ss = results_bss.x + + # Generate other steady-state values and Euler equations + K_ss, K_cnstr = get_K(b_ss) + r_params = (A, alpha, delta) + r_ss = get_r(r_params, K_ss, L) + w_params = (A, alpha) + w_ss = get_w(w_params, K_ss, L) + c_ss, c_cnstr = get_cvec(r_ss, w_ss, b_ss, nvec) + Y_params = (A, alpha) + Y_ss = get_Y(Y_params, K_ss, L) + C_ss = get_C(c_ss) + b_err_params = (beta, sigma) + EulErr_ss = get_b_errors( + b_err_params, r_ss, c_ss, c_cnstr, EulDiff) + RCerr_ss = Y_ss - C_ss - delta * K_ss + + ss_time = time.time() - start_time + + ss_output = { + 'b_ss': b_ss, 'c_ss': c_ss, 'w_ss': w_ss, 'r_ss': r_ss, + 'K_ss': K_ss, 'Y_ss': Y_ss, 'C_ss': C_ss, + 'EulErr_ss': EulErr_ss, 'RCerr_ss': RCerr_ss, + 'ss_time': ss_time} + print('b_ss is: ', b_ss) + print('c_ss is: ', c_ss) + print('Euler errors are: ', EulErr_ss) + print('Resource constraint error is: ', RCerr_ss) + + # Print SS computation time + print_time(ss_time, 'SS') + + if graphs: + ''' + ---------------------------------------------------------------- + cur_path = string, path name of current directory + output_fldr = string, folder in current path to save files + output_dir = string, total path of images folder + output_path = string, path of file name of figure to be saved + S = integer >= 3, number of periods in a life + age_pers = (S,) vector, ages from 1 to S + ---------------------------------------------------------------- + ''' + # Create directory if images directory does not already exist + cur_path = os.path.split(os.path.abspath(__file__))[0] + output_fldr = "images" + output_dir = os.path.join(cur_path, output_fldr) + if not os.access(output_dir, os.F_OK): + os.makedirs(output_dir) + + # Plot steady-state consumption and savings distributions + S = nvec.shape[0] + age_pers = np.arange(1, S + 1) + fig, ax = plt.subplots() + plt.plot(age_pers, c_ss, marker='D', label='Consumption') + plt.plot(age_pers, np.hstack((0, b_ss)), marker='D', + label='Savings') + # for the minor ticks, use no labels; default NullFormatter + minorLocator = MultipleLocator(1) + ax.xaxis.set_minor_locator(minorLocator) + plt.grid(visible=True, which='major', color='0.65', linestyle='-') + plt.title('Steady-state consumption and savings', fontsize=20) + plt.xlabel(r'Age $s$') + plt.ylabel(r'Units of consumption') + plt.xlim((0.8, S + 0.2)) + plt.ylim((-0.02, 1.15 * (c_ss.max()))) + plt.legend(loc='center right') + output_path = os.path.join(output_dir, "SS_bc") + plt.savefig(output_path) + # plt.show() + + return ss_output diff --git a/docs/book/content/UNtutorial/solutions/3perOG/TPI.py b/docs/book/content/UNtutorial/solutions/3perOG/TPI.py new file mode 100644 index 0000000..7f53bc3 --- /dev/null +++ b/docs/book/content/UNtutorial/solutions/3perOG/TPI.py @@ -0,0 +1,640 @@ +''' +------------------------------------------------------------------------ +This module contains the functions used to solve the time path iteration +non-steady-state equilibrium for the model with 3-period lived agents +and exogenous labor from Chapter 5 of the OG textbook. + +This Python module calls the following function(s): + get_path() + get_cvec_lf() + LfEulerSys() + paths_life() + get_cbepath() + get_TPI() + ss.print_time() + + get_cvec() + get_L() + get_K() + get_w() + get_r() + get_Y() + get_C() + feasible() + EulerSys() + get_b_errors() + get_SS() +------------------------------------------------------------------------ +''' +# Import Packages +import time +import numpy as np +import scipy.optimize as opt +import SS as ss +import matplotlib +import matplotlib.pyplot as plt +from matplotlib.ticker import MultipleLocator +from mpl_toolkits.mplot3d import Axes3D +import sys +import os + +''' +------------------------------------------------------------------------ + Functions +------------------------------------------------------------------------ +''' + + +def get_path(x1, xT, T, spec): + ''' + -------------------------------------------------------------------- + This function generates a path from point x1 to point xT such that + that the path x is a linear or quadratic function of time t. + + linear: x = d*t + e + quadratic: x = a*t^2 + b*t + c + + The identifying assumptions for quadratic are the following: + + (1) x1 is the value at time t=0: x1 = c + (2) xT is the value at time t=T-1: xT = a*(T-1)^2 + b*(T-1) + c + (3) the slope of the path at t=T-1 is 0: 0 = 2*a*(T-1) + b + -------------------------------------------------------------------- + INPUTS: + x1 = scalar, initial value of the function x(t) at t=0 + xT = scalar, value of the function x(t) at t=T-1 + T = integer >= 3, number of periods of the path + spec = string, "linear" or "quadratic" + + OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION: None + + OBJECTS CREATED WITHIN FUNCTION: + cc = scalar, constant coefficient in quadratic function + bb = scalar, coefficient on t in quadratic function + aa = scalar, coefficient on t^2 in quadratic function + xpath = (T,) vector, parabolic xpath from x1 to xT + + FILES CREATED BY THIS FUNCTION: None + + RETURNS: xpath + -------------------------------------------------------------------- + ''' + if spec == "linear": + xpath = np.linspace(x1, xT, T) + elif spec == "quadratic": + cc = x1 + bb = 2 * (xT - x1) / (T - 1) + aa = (x1 - xT) / ((T - 1) ** 2) + xpath = aa * (np.arange(0, T) ** 2) + (bb * np.arange(0, T)) + cc + + return xpath + + +def get_cvec_lf(rpath, wpath, nvec, bvec): + ''' + -------------------------------------------------------------------- + Generates vector of remaining lifetime consumptions from individual + savings, and the time path of interest rates and the real wages, + where p is an integer in [2, S] representing the remaining periods + of life + -------------------------------------------------------------------- + INPUTS: + rpath = (p,) vector, remaining interest rates + wpath = (p,) vector, remaining wages + nvec = (p,) vector, remaining exogenous labor supply + bvec = (p,) vector, remaining savings including initial savings + + OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION: None + + OBJECTS CREATED WITHIN FUNCTION: + b_s = (p,) vector, bvec + b_sp1 = (p,) vector, last p-1 elements of bvec and 0 in last + element + cvec = (p,) vector, remaining consumption by age c_s + c_cnstr = (p,) Boolean vector, =True if element c_s <= 0 + + FILES CREATED BY THIS FUNCTION: None + + RETURNS: cvec, c_cnstr + -------------------------------------------------------------------- + ''' + b_s = bvec + b_sp1 = np.append(bvec[1:], [0]) + cvec = (1 + rpath) * b_s + wpath * nvec - b_sp1 + if cvec.min() <= 0: + print('get_cvec_lf() warning: distribution of savings and/or ' + + 'parameters created c<=0 for some agent(s)') + c_cnstr = cvec <= 0 + return cvec, c_cnstr + + +def LfEulerSys(bvec, *args): + ''' + -------------------------------------------------------------------- + Generates vector of all Euler errors for a given bvec, which errors + characterize all optimal lifetime decisions, where p is an integer + in [2, S] representing the remaining periods of life + -------------------------------------------------------------------- + INPUTS: + bvec = (p-1,) vector, remaining lifetime savings decisions + where p is the number of remaining periods + args = length 7 tuple, (beta, sigma, beg_wealth, nvec, rpath, + wpath, EulDiff) + beta = scalar in [0,1), discount factor + sigma = scalar > 0, coefficient of relative risk aversion + beg_wealth = scalar, wealth at the beginning of first age + nvec = (p,) vector, remaining exogenous labor supply + rpath = (p,) vector, interest rates over remaining life + wpath = (p,) vector, wages rates over remaining life + + OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION: + get_cvec_lf() + c5ssf.get_b_errors() + + OBJECTS CREATED WITHIN FUNCTION: + bvec2 = (p,) vector, remaining savings including initial + savings + cvec = (p,) vector, remaining lifetime consumption + levels implied by bvec2 + c_cnstr = (p,) Boolean vector, =True if c_{s,t}<=0 + b_err_params = length 2 tuple, (beta, sigma) + b_err_vec = (p-1,) vector, Euler errors from lifetime + consumption vector + + FILES CREATED BY THIS FUNCTION: None + + RETURNS: b_err_vec + -------------------------------------------------------------------- + ''' + beta, sigma, beg_wealth, nvec, rpath, wpath, EulDiff = args + bvec2 = np.append(beg_wealth, bvec) + cvec, c_cnstr = get_cvec_lf(rpath, wpath, nvec, bvec2) + b_err_params = (beta, sigma) + b_err_vec = ss.get_b_errors(b_err_params, rpath[1:], cvec, + c_cnstr, EulDiff) + return b_err_vec + + +def paths_life(params, beg_age, beg_wealth, nvec, rpath, wpath, + b_init): + ''' + -------------------------------------------------------------------- + Solve for the remaining lifetime savings decisions of an individual + who enters the model at age beg_age, with corresponding initial + wealth beg_wealth. Variable p is an integer in [2, S] representing + the remaining periods of life. + -------------------------------------------------------------------- + INPUTS: + params = length 5 tuple, (S, beta, sigma, TPI_tol, EulDiff) + S = integer in [3,80], number of periods an individual + lives + beta = scalar in (0,1), discount factor for each model + period + sigma = scalar > 0, coefficient of relative risk aversion + TPI_tol = scalar > 0, tolerance level for fsolve's in TPI + EulDiff = Boolean, =True if want difference version of Euler + errors beta*(1+r)*u'(c2) - u'(c1), =False if want + ratio version [beta*(1+r)*u'(c2)]/[u'(c1)] - 1 + beg_age = integer in [1,S-1], beginning age of remaining life + beg_wealth = scalar, beginning wealth at beginning age + nvec = (p,) vector, remaining exogenous labor supplies + rpath = (p,) vector, remaining lifetime interest rates + wpath = (p,) vector, remaining lifetime wages + b_init = (p-1,) vector, initial guess for remaining lifetime + savings + + OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION: + LfEulerSys() + get_cvec_lf() + c4ssf.get_b_errors() + + OBJECTS CREATED WITHIN FUNCTION: + p = integer in [2,S], remaining periods in life + b_guess = (p-1,) vector, initial guess for lifetime savings + decisions + eullf_objs = length 7 tuple, (beta, sigma, beg_wealth, nvec, + rpath, wpath, EulDiff) + bpath = (p-1,) vector, optimal remaining lifetime savings + decisions + cpath = (p,) vector, optimal remaining lifetime + consumption decisions + c_cnstr = (p,) boolean vector, =True if c_p <= 0, + b_err_params = length 2 tuple, (beta, sigma) + b_err_vec = (p-1,) vector, Euler errors associated with + optimal savings decisions + + FILES CREATED BY THIS FUNCTION: None + + RETURNS: bpath, cpath, b_err_vec + -------------------------------------------------------------------- + ''' + S, beta, sigma, TPI_tol, EulDiff = params + p = int(S - beg_age + 1) + if beg_age == 1 and beg_wealth != 0: + sys.exit("Beginning wealth is nonzero for age s=1.") + if len(rpath) != p: + sys.exit("Beginning age and length of rpath do not match.") + if len(wpath) != p: + sys.exit("Beginning age and length of wpath do not match.") + if len(nvec) != p: + sys.exit("Beginning age and length of nvec do not match.") + b_guess = 1.01 * b_init + eullf_objs = (beta, sigma, beg_wealth, nvec, rpath, wpath, EulDiff) + bpath = opt.fsolve(LfEulerSys, b_guess, args=(eullf_objs), + xtol=TPI_tol) + cpath, c_cnstr = get_cvec_lf(rpath, wpath, nvec, + np.append(beg_wealth, bpath)) + b_err_params = (beta, sigma) + b_err_vec = ss.get_b_errors(b_err_params, rpath[1:], cpath, + c_cnstr, EulDiff) + return bpath, cpath, b_err_vec + + +def get_cbepath(params, rpath, wpath): + ''' + -------------------------------------------------------------------- + Generates matrices for the time path of the distribution of + individual savings, individual consumption, and the Euler errors + associated with the savings decisions. + -------------------------------------------------------------------- + INPUTS: + params = length 9 tuple, + (S, T, beta, sigma, nvec, bvec1, b_ss, TPI_tol, EulDiff) + S = integer in [3,80], number of periods an individual lives + T = integer > S, number of time periods until steady state + beta = scalar in (0,1), discount factor for each model period + sigma = scalar > 0, coefficient of relative risk aversion + nvec = (S,) vector, exogenous labor supply n_s + bvec1 = (S-1,) vector, initial period savings distribution + b_ss = (S-1,) vector, steady-state savings distribution + TPI_tol = scalar > 0, tolerance level for fsolve's in TPI + EulDiff = Boolean, =True if want difference version of Euler errors + beta*(1+r)*u'(c2) - u'(c1), =False if want ratio version + [beta*(1+r)*u'(c2)]/[u'(c1)] - 1 + rpath = (T+S-2,) vector, equilibrium time path of interest rate + wpath = (T+S-2,) vector, equilibrium time path of the real wage + + OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION: + paths_life() + + OBJECTS CREATED WITHIN FUNCTION: + cpath = (S, T+S-2) matrix, time path of the distribution of + consumption + bpath = (S-1, T+S-2) matrix, time path of the distribution of + savings + EulErrPath = (S-1, T+S-2) matrix, time path of Euler errors + pl_params = length 5 tuple, (S, beta, sigma, TPI_tol, EulDiff) + p = integer >= 2, index representing number of periods + remaining in a lifetime, used to solve incomplete + lifetimes + b_guess = (p-1,) vector, initial guess for remaining lifetime + savings, taken from previous cohort's choices + bveclf = (p-1,) vector, optimal remaining lifetime savings + decisions + cveclf = (p,) vector, optimal remaining lifetime consumption + decisions + b_err_veclf = (p-1,) vector, Euler errors associated with + optimal remaining lifetime savings decisions + DiagMaskb = (p-1, p-1) Boolean identity matrix + DiagMaskc = (p, p) Boolean identity matrix + + FILES CREATED BY THIS FUNCTION: None + + RETURNS: cpath, bpath, EulErrPath + -------------------------------------------------------------------- + ''' + S, T, beta, sigma, nvec, bvec1, b_ss, TPI_tol, EulDiff = params + cpath = np.zeros((S, T + S - 2)) + bpath = np.append(bvec1.reshape((S - 1, 1)), + np.zeros((S - 1, T + S - 3)), axis=1) + EulErrPath = np.zeros((S - 1, T + S - 2)) + # Solve the incomplete remaining lifetime decisions of agents alive + # in period t=1 but not born in period t=1 + cpath[S - 1, 0] = ((1 + rpath[0]) * bvec1[S - 2] + + wpath[0] * nvec[S - 1]) + pl_params = (S, beta, sigma, TPI_tol, EulDiff) + for p in range(2, S): + b_guess = np.diagonal(bpath[S - p:, :p - 1]) + bveclf, cveclf, b_err_veclf = paths_life( + pl_params, S - p + 1, bvec1[S - p - 1], nvec[-p:], + rpath[:p], wpath[:p], b_guess) + # Insert the vector lifetime solutions diagonally (twist donut) + # into the cpath, bpath, and EulErrPath matrices + DiagMaskb = np.eye(p - 1, dtype=bool) + DiagMaskc = np.eye(p, dtype=bool) + bpath[S - p:, 1:p] = DiagMaskb * bveclf + bpath[S - p:, 1:p] + cpath[S - p:, :p] = DiagMaskc * cveclf + cpath[S - p:, :p] + EulErrPath[S - p:, 1:p] = (DiagMaskb * b_err_veclf + + EulErrPath[S - p:, 1:p]) + # Solve for complete lifetime decisions of agents born in periods + # 1 to T and insert the vector lifetime solutions diagonally (twist + # donut) into the cpath, bpath, and EulErrPath matrices + DiagMaskb = np.eye(S - 1, dtype=bool) + DiagMaskc = np.eye(S, dtype=bool) + for t in range(1, T): # Go from periods 1 to T-1 + b_guess = np.diagonal(bpath[:, t - 1:t + S - 2]) + bveclf, cveclf, b_err_veclf = paths_life( + pl_params, 1, 0, nvec, rpath[t - 1:t + S - 1], + wpath[t - 1:t + S - 1], b_guess) + # Insert the vector lifetime solutions diagonally (twist donut) + # into the cpath, bpath, and EulErrPath matrices + bpath[:, t:t + S - 1] = (DiagMaskb * bveclf + + bpath[:, t:t + S - 1]) + cpath[:, t - 1:t + S - 1] = (DiagMaskc * cveclf + + cpath[:, t - 1:t + S - 1]) + EulErrPath[:, t:t + S - 1] = (DiagMaskb * b_err_veclf + + EulErrPath[:, t:t + S - 1]) + + return cpath, bpath, EulErrPath + + +def get_TPI(params, bvec1, graphs): + ''' + -------------------------------------------------------------------- + Generates steady-state time path for all endogenous objects from + initial state (K1, Gamma1) to the steady state. + -------------------------------------------------------------------- + INPUTS: + params = length 17 tuple, (S, T, beta, sigma, nvec, L, A, + alpha, delta, b_ss, K_ss, C_ss, maxiter_TPI, + mindist_TPI, xi, TPI_tol, EulDiff) + S = integer in [3,80], number of periods an individual + lives + T = integer > S, number of time periods until steady + state + beta = scalar in (0,1), discount factor for model period + sigma = scalar > 0, coefficient of relative risk aversion + nvec = (S,) vector, exogenous labor supply n_{s,t} + L = scalar > 0, exogenous aggregate labor + A = scalar > 0, total factor productivity parameter in + firms' production function + alpha = scalar in (0,1), capital share of income + delta = scalar in [0,1], model-period depreciation rate of + capital + b_ss = (S-1,) vector, steady-state distribution of savings + K_ss = scalar > 0, steady-state aggregate capital stock + C_ss = scalar > 0, steady-state aggregate consumption + maxiter_TPI = integer >= 1, Maximum number of iterations for TPI + mindist_TPI = scalar > 0, Convergence criterion for TPI + xi = scalar in (0,1], TPI path updating parameter + TPI_tol = scalar > 0, tolerance level for fsolve's in TPI + EulDiff = Boolean, =True if want difference version of Euler + errors beta*(1+r)*u'(c2) - u'(c1), =False if want + ratio version [beta*(1+r)*u'(c2)]/[u'(c1)] - 1 + bvec1 = (S-1,) vector, initial period savings distribution + graphs = Boolean, =True if want graphs of TPI objects + + OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION: + c5ssf.get_K() + get_path() + c5ssf.get_r() + c5ssf.get_w() + get_cbepath() + c5ssf.get_Y() + c5ssf.get_C() + c5ssf.print_time() + + + OBJECTS CREATED WITHIN FUNCTION: + start_time = scalar, current processor time in seconds (float) + K1 = scalar > 0, initial aggregate capital stock + K1_cnstr = Boolean, =True if K1 <= 0 + Kpath_init = (T+S-2,) vector, initial guess for the time path + of the aggregate capital stock + Lpath = (T+S-2,) vector, exogenous time path for aggregate + labor + iter_TPI = integer >= 0, current iteration of TPI + dist_TPI = scalar >= 0, distance measure for fixed point + Kpath_new = (T+S-2,) vector, new path of the aggregate capital + stock implied by household and firm optimization + r_params = length 3 tuple, (A, alpha, delta) + w_params = length 2 tuple, (A, alpha) + cbe_params = length 9 tuple, (S, T, beta, sigma, nvec, bvec1, b_ss, + TPI_tol, EulDiff) + rpath = (T+S-2,) vector, time path of the interest rate + wpath = (T+S-2,) vector, time path of the wage + cpath = (S, T+S-2) matrix, time path values of distribution of + consumption c_{s,t} + bpath = (S-1, T+S-2) matrix, time path of distribution of + individual savings b_{s,t} + EulErrPath = (S-1, T+S-2) matrix, time path of individual Euler + errors corresponding to individual savings b_{s,t} + (first column is zeros) + Kpath_cnstr = (T+S-2,) Boolean vector, =True if K_t<=0 + Kpath = (T+S-2,) vector, equilibrium time path of aggregate + capital stock K_t + Y_params = length 2 tuple, (A, alpha) + Ypath = (T+S-2,) vector, equilibrium time path of aggregate + output (GDP) Y_t + Cpath = (T+S-2,) vector, equilibrium time path of aggregate + consumption C_t + RCerrPath = (T+S-1,) vector, equilibrium time path of the resource + constraint error: Y_t - C_t - K_{t+1} + (1-delta)*K_t + tpi_time = scalar, time to compute TPI solution (seconds) + tpi_output = length 10 dictionary, {bpath, cpath, wpath, rpath, + Kpath, Ypath, Cpath, EulErrPath, RCerrPath, tpi_time} + + FILES CREATED BY THIS FUNCTION: + Kpath.png + Ypath.png + Cpath.png + wpath.png + rpath.png + bpath.png + cpath.png + + RETURNS: tpi_output + -------------------------------------------------------------------- + ''' + start_time = time.time() + (S, T, beta, sigma, nvec, L, A, alpha, delta, b_ss, K_ss, C_ss, + maxiter_TPI, mindist_TPI, xi, TPI_tol, EulDiff) = params + K1, K1_cnstr = ss.get_K(bvec1) + + # Create time paths for K and L + Kpath_init = np.zeros(T + S - 2) + Kpath_init[:T] = get_path(K1, K_ss, T, "quadratic") + Kpath_init[T:] = K_ss + Lpath = L * np.ones(T + S - 2) + + iter_TPI = int(0) + dist_TPI = 10. + Kpath_new = Kpath_init.copy() + r_params = (A, alpha, delta) + w_params = (A, alpha) + cbe_params = (S, T, beta, sigma, nvec, bvec1, b_ss, TPI_tol, + EulDiff) + + while (iter_TPI < maxiter_TPI) and (dist_TPI >= mindist_TPI): + iter_TPI += 1 + Kpath_init = xi * Kpath_new + (1 - xi) * Kpath_init + rpath = ss.get_r(r_params, Kpath_init, Lpath) + wpath = ss.get_w(w_params, Kpath_init, Lpath) + cpath, bpath, EulErrPath = get_cbepath(cbe_params, rpath, wpath) + Kpath_new = np.zeros(T + S - 2) + Kpath_new[:T], Kpath_cnstr = ss.get_K(bpath[:, :T]) + Kpath_new[T:] = K_ss * np.ones(S - 2) + Kpath_cnstr = np.append(Kpath_cnstr, + np.zeros(S - 2, dtype=bool)) + Kpath_new[Kpath_cnstr] = 0.1 + # Check the distance of Kpath_new1 + dist_TPI = ((Kpath_new[1:T] - Kpath_init[1:T]) ** 2).sum() + # dist_TPI = np.absolute((Kpath_new[1:T] - Kpath_init[1:T]) / + # Kpath_init[1:T]).max() + print('iter: ', iter_TPI, ', dist: ', dist_TPI, + ',max Eul err: ', np.absolute(EulErrPath).max()) + + if iter_TPI == maxiter_TPI and dist_TPI > mindist_TPI: + print('TPI reached maxiter and did not converge.') + elif iter_TPI == maxiter_TPI and dist_TPI <= mindist_TPI: + print('TPI converged in the last iteration. ' + + 'Should probably increase maxiter_TPI.') + Kpath = Kpath_new + Y_params = (A, alpha) + Ypath = ss.get_Y(Y_params, Kpath, Lpath) + Cpath = np.zeros(T + S - 2) + Cpath[:T - 1] = ss.get_C(cpath[:, :T - 1]) + Cpath[T - 1:] = C_ss * np.ones(S - 1) + RCerrPath = (Ypath[:-1] - Cpath[:-1] - Kpath[1:] + + (1 - delta) * Kpath[:-1]) + tpi_time = time.time() - start_time + + tpi_output = { + 'bpath': bpath, 'cpath': cpath, 'wpath': wpath, 'rpath': rpath, + 'Kpath': Kpath, 'Ypath': Ypath, 'Cpath': Cpath, + 'EulErrPath': EulErrPath, 'RCerrPath': RCerrPath, + 'tpi_time': tpi_time} + + # Print TPI computation time + ss.print_time(tpi_time, 'TPI') + + if graphs: + ''' + ---------------------------------------------------------------- + cur_path = string, path name of current directory + output_fldr = string, folder in current path to save files + output_dir = string, total path of images folder + output_path = string, path of file name of figure to be saved + tvec = (T+S-2,) vector, time period vector + tgridTm1 = (T-1,) vector, time period vector to T-1 + tgridT = (T,) vector, time period vector to T-1 + sgrid = (S,) vector, all ages from 1 to S + sgrid2 = (S-1,) vector, all ages from 2 to S + tmatb = (2, 18) matrix, time periods for all savings + decisions ages (S-1) and time periods (T) + smatb = (2, 18) matrix, ages for all savings decision ages + (S-1) and time periods (T) + tmatc = (3, 17) matrix, time periods for all consumption + decisions ages (S) and time periods (T-1) + smatc = (3, 17) matrix, ages for all consumption decisions + ages (S) and time periods (T-1) + ---------------------------------------------------------------- + ''' + # Create directory if images directory does not already exist + cur_path = os.path.split(os.path.abspath(__file__))[0] + output_fldr = "images" + output_dir = os.path.join(cur_path, output_fldr) + if not os.access(output_dir, os.F_OK): + os.makedirs(output_dir) + + # Plot time path of aggregate capital stock + tvec = np.linspace(1, T + S - 2, T + S - 2) + minorLocator = MultipleLocator(1) + fig, ax = plt.subplots() + plt.plot(tvec, Kpath, marker='D') + # for the minor ticks, use no labels; default NullFormatter + ax.xaxis.set_minor_locator(minorLocator) + plt.grid(visible=True, which='major', color='0.65', linestyle='-') + plt.title('Time path for aggregate capital stock K') + plt.xlabel(r'Period $t$') + plt.ylabel(r'Aggregate capital $K_{t}$') + output_path = os.path.join(output_dir, "Kpath") + plt.savefig(output_path) + # plt.show() + + # Plot time path of aggregate output (GDP) + fig, ax = plt.subplots() + plt.plot(tvec, Ypath, marker='D') + # for the minor ticks, use no labels; default NullFormatter + ax.xaxis.set_minor_locator(minorLocator) + plt.grid(visible=True, which='major', color='0.65', linestyle='-') + plt.title('Time path for aggregate output (GDP) Y') + plt.xlabel(r'Period $t$') + plt.ylabel(r'Aggregate output $Y_{t}$') + output_path = os.path.join(output_dir, "Ypath") + plt.savefig(output_path) + # plt.show() + + # Plot time path of aggregate consumption + fig, ax = plt.subplots() + plt.plot(tvec, Cpath, marker='D') + # for the minor ticks, use no labels; default NullFormatter + ax.xaxis.set_minor_locator(minorLocator) + plt.grid(visible=True, which='major', color='0.65', linestyle='-') + plt.title('Time path for aggregate consumption C') + plt.xlabel(r'Period $t$') + plt.ylabel(r'Aggregate consumption $C_{t}$') + output_path = os.path.join(output_dir, "C_aggr_path") + plt.savefig(output_path) + # plt.show() + + # Plot time path of real wage + fig, ax = plt.subplots() + plt.plot(tvec, wpath, marker='D') + # for the minor ticks, use no labels; default NullFormatter + ax.xaxis.set_minor_locator(minorLocator) + plt.grid(visible=True, which='major', color='0.65', linestyle='-') + plt.title('Time path for real wage w') + plt.xlabel(r'Period $t$') + plt.ylabel(r'Real wage $w_{t}$') + output_path = os.path.join(output_dir, "wpath") + plt.savefig(output_path) + # plt.show() + + # Plot time path of real interest rate + fig, ax = plt.subplots() + plt.plot(tvec, rpath, marker='D') + # for the minor ticks, use no labels; default NullFormatter + ax.xaxis.set_minor_locator(minorLocator) + plt.grid(visible=True, which='major', color='0.65', linestyle='-') + plt.title('Time path for real interest rate r') + plt.xlabel(r'Period $t$') + plt.ylabel(r'Real interest rate $r_{t}$') + output_path = os.path.join(output_dir, "rpath") + plt.savefig(output_path) + # plt.show() + + # Plot time path of individual savings distribution + tgridT = np.linspace(1, T, T) + sgrid2 = np.linspace(2, S, S - 1) + tmatb, smatb = np.meshgrid(tgridT, sgrid2) + cmap_bp = matplotlib.cm.get_cmap('summer') + fig, ax = plt.subplots(subplot_kw={"projection": "3d"}) + ax.set_xlabel(r'period-$t$') + ax.set_ylabel(r'age-$s$') + ax.set_zlabel(r'individual savings $b_{s,t}$') + strideval = max(int(1), int(round(S / 10))) + ax.plot_surface(tmatb, smatb, bpath[:, :T], rstride=strideval, + cstride=strideval, cmap=cmap_bp) + output_path = os.path.join(output_dir, "bpath") + plt.savefig(output_path) + # plt.show() + + # Plot time path of individual consumption distribution + tgridTm1 = np.linspace(1, T - 1, T - 1) + sgrid = np.linspace(1, S, S) + tmatc, smatc = np.meshgrid(tgridTm1, sgrid) + cmap_cp = matplotlib.cm.get_cmap('summer') + fig, ax = plt.subplots(subplot_kw={"projection": "3d"}) + ax.set_xlabel(r'period-$t$') + ax.set_ylabel(r'age-$s$') + ax.set_zlabel(r'individual consumption $c_{s,t}$') + strideval = max(int(1), int(round(S / 10))) + ax.plot_surface(tmatc, smatc, cpath[:, :T - 1], + rstride=strideval, cstride=strideval, + cmap=cmap_cp) + output_path = os.path.join(output_dir, "cpath") + plt.savefig(output_path) + # plt.show() + + return tpi_output diff --git a/docs/book/content/UNtutorial/solutions/3perOG/execute.py b/docs/book/content/UNtutorial/solutions/3perOG/execute.py new file mode 100755 index 0000000..41afb43 --- /dev/null +++ b/docs/book/content/UNtutorial/solutions/3perOG/execute.py @@ -0,0 +1,200 @@ +''' +------------------------------------------------------------------------ +This program runs the steady state solver as well as the time path +iteration solution for the model with 3-period lived agents and +exogenous labor from Chapter 5 of the OG textbook. + +This Python script imports the following module(s): + SS.py + TPI.py + +This Python script calls the following function(s): + ss.feasible() + ss.get_SS() + ss.get_K + tpi.get_TPI() +------------------------------------------------------------------------ +''' +# Import packages +import numpy as np +import SS as ss +import TPI as tpi + +''' +------------------------------------------------------------------------ +Declare parameters +------------------------------------------------------------------------ +S = integer in [3,80], number of periods an individual lives +beta_annual = scalar in (0,1), discount factor for one year +beta = scalar in (0,1), discount factor for each model period +sigma = scalar > 0, coefficient of relative risk aversion +nvec = [S,] vector, exogenous labor supply n_{s,t} +L = scalar > 0, exogenous aggregate labor +A = scalar > 0, total factor productivity parameter in firms' + production function +alpha = scalar in (0,1), capital share of income +delta_annual = scalar in [0,1], one-year depreciation rate of capital +delta = scalar in [0,1], model-period depreciation rate of + capital +SS_tol = scalar > 0, tolerance level for steady-state fsolve +SS_graphs = boolean, =True if want graphs of steady-state objects +T = integer > S, number of time periods until steady state +TPI_solve = boolean, =True if want to solve TPI after solving SS +TPI_tol = scalar > 0, tolerance level for fsolve's in TPI +maxiter_TPI = integer >= 1, Maximum number of iterations for TPI +mindist_TPI = scalar > 0, Convergence criterion for TPI +xi = scalar in (0,1], TPI path updating parameter +TPI_graphs = Boolean, =True if want graphs of TPI objects +EulDiff = Boolean, =True if want difference version of Euler errors + beta*(1+r)*u'(c2) - u'(c1), =False if want ratio version + [beta*(1+r)*u'(c2)]/[u'(c1)] - 1 +------------------------------------------------------------------------ +''' +# Household parameters +S = int(3) +beta_annual = .96 +beta = beta_annual ** 20 +sigma = 3.0 +nvec = np.array([1.0, 1.0, 0.2]) +L = nvec.sum() +# Firm parameters +A = 1.0 +alpha = .35 +delta_annual = .05 +delta = 1 - ((1 - delta_annual) ** 20) +# SS parameters +SS_tol = 1e-13 +SS_graphs = True +# TPI parameters +T = int(round(6 * S)) +TPI_solve = True +TPI_tol = 1e-13 +maxiter_TPI = 200 +mindist_TPI = 1e-13 +xi = 0.99 +TPI_graphs = True +# Overall parameters +EulDiff = False + +''' +------------------------------------------------------------------------ +Check feasibility +------------------------------------------------------------------------ +f_params = length 4 tuple, (nvec, A, alpha, delta) +bvec_guess1 = (2,) vector, guess for steady-state bvec (b1, b2) +b_cnstr = (2,) Boolean vector, =True if b_s causes negative + consumption c_s <= 0 or negative aggregate capital stock + K <= 0 +c_cnstr = (3,) Boolean vector, =True for elements of negative + consumption c_s <= 0 +K_cnstr = Boolean, =True if K <= 0 +bvec_guess2 = (2,) vector, guess for steady-state bvec (b1, b2) +bvec_guess3 = (2,) vector, guess for steady-state bvec (b1, b2) + +------------------------------------------------------------------------ +''' +f_params = (nvec, A, alpha, delta) + +bvec_guess1 = np.array([1.0, 1.2]) +b_cnstr, c_cnstr, K_cnstr = ss.feasible(f_params, bvec_guess1) +print('bvec_guess1', bvec_guess1) +print('c_cnstr', c_cnstr) +print('K_cnstr', K_cnstr) + +bvec_guess2 = np.array([0.06, -0.001]) +b_cnstr, c_cnstr, K_cnstr = ss.feasible(f_params, bvec_guess2) +print('bvec_guess2', bvec_guess2) +print('c_cnstr', c_cnstr) +print('K_cnstr', K_cnstr) + +bvec_guess3 = np.array([0.1, 0.1]) +b_cnstr, c_cnstr, K_cnstr = ss.feasible(f_params, bvec_guess3) +print('bvec_guess3', bvec_guess3) +print('c_cnstr', c_cnstr) +print('K_cnstr', K_cnstr) + +''' +------------------------------------------------------------------------ +Run the steady-state solution +------------------------------------------------------------------------ +bvec_guess = (2,) vector, initial guess for steady-state bvec (b1, b2) +f_params = length 4 tuple, (nvec, A, alpha, delta) +b_cnstr = (2,) Boolean vector, =True if b_s causes negative + consumption c_s <= 0 or negative aggregate capital stock + K <= 0 +c_cnstr = (3,) Boolean vector, =True for elements of negative + consumption c_s <= 0 +K_cnstr = Boolean, =True if K <= 0 +ss_params = length 9 tuple, + (beta, sigma, nvec, L, A, alpha, delta, SS_tol, EulDiff) +ss_output = length 10 dictionary, {b_ss, c_ss, w_ss, r_ss, K_ss, Y_ss, + C_ss, EulErr_ss, RCerr_ss, ss_time} +------------------------------------------------------------------------ +''' +print('BEGIN EQUILIBRIUM STEADY-STATE COMPUTATION') +bvec_guess = np.array([0.1, 0.1]) +f_params = (nvec, A, alpha, delta) +b_cnstr, c_cnstr, K_cnstr = ss.feasible(f_params, bvec_guess) +if not K_cnstr and not c_cnstr.max(): + ss_params = (beta, sigma, nvec, L, A, alpha, delta, SS_tol, + EulDiff) + ss_output = ss.get_SS(ss_params, bvec_guess, SS_graphs) + + beta2 = 0.55 + ss_params2 = (beta2, sigma, nvec, L, A, alpha, delta, SS_tol, + EulDiff) + ss_output2 = ss.get_SS(ss_params2, bvec_guess, False) + print('c2_ss = ', ss_output2['c_ss'], ', diff is ', + ss_output2['c_ss'] - ss_output['c_ss']) + print('b2_ss = ', ss_output2['b_ss'], ', diff is ', + ss_output2['b_ss'] - ss_output['b_ss']) + print('w2_ss = ', ss_output2['w_ss'], ', diff is ', + ss_output2['w_ss'] - ss_output['w_ss']) + print('r2_ss = ', ss_output2['r_ss'], ', diff is ', + ss_output2['r_ss'] - ss_output['r_ss']) + print('K2_ss = ', ss_output2['K_ss'], ', diff is ', + ss_output2['K_ss'] - ss_output['K_ss']) + print('Y2_ss = ', ss_output2['Y_ss'], ', diff is ', + ss_output2['Y_ss'] - ss_output['Y_ss']) + print('C2_ss = ', ss_output2['C_ss'], ', diff is ', + ss_output2['C_ss'] - ss_output['C_ss']) + print('EulErr2_ss = ', ss_output2['EulErr_ss'], ', diff is ', + ss_output2['EulErr_ss'] - ss_output['EulErr_ss']) + print('RCerr2_ss = ', ss_output2['RCerr_ss'], ', diff is ', + ss_output2['RCerr_ss'] - ss_output['RCerr_ss']) +else: + print("Initial guess for SS bvec does not satisfy K>0 or c_s>0.") + +''' +------------------------------------------------------------------------ +Run the time path iteration (TPI) solution +------------------------------------------------------------------------ +b_ss = (2,) vector, steady-state savings distribution +K_ss = scalar > 0, steady-state aggregate capital stock +C_ss = scalar > 0, steady-state aggregate consumption +bvec1 = (2,) vector, initial period savings distribution +K1 = scalar, initial period aggregate capital stock +K_constr_tp1 = Boolean, =True if K1 <= 0 +tpi_params = length 17 tuple, (S, T, beta, sigma, nvec, L, A, alpha, + delta, b_ss, K_ss, C_ss, maxiter_TPI, mindist_TPI, xi, + TPI_tol, EulDiff) +tpi_output = length 10 dictionary, {bpath, cpath, wpath, rpath, Kpath, + Ypath, Cpath, EulErrPath, RCerrPath, tpi_time} +------------------------------------------------------------------------ +''' +if TPI_solve: + print('BEGIN EQUILIBRIUM TIME PATH COMPUTATION') + b_ss = ss_output['b_ss'] + K_ss = ss_output['K_ss'] + C_ss = ss_output['C_ss'] + bvec1 = np.array([0.8 * b_ss[0], 1.1 * b_ss[1]]) + # Make sure init. period distribution is feasible in terms of K + K1, K_constr_tpi1 = ss.get_K(bvec1) + if K_constr_tpi1: + print('Initial savings distribution is not feasible because ' + + 'K1<=0. Some element(s) of bvec1 must increase.') + else: + tpi_params = (S, T, beta, sigma, nvec, L, A, alpha, delta, b_ss, + K_ss, C_ss, maxiter_TPI, mindist_TPI, xi, TPI_tol, + EulDiff) + tpi_output = tpi.get_TPI(tpi_params, bvec1, TPI_graphs) From 4f2cd7677d5cc1a60812d7e6a20393d210d2f39a Mon Sep 17 00:00:00 2001 From: Richard Evans Date: Wed, 23 Oct 2024 00:22:25 +0800 Subject: [PATCH 2/4] Fixed link type --- docs/book/content/UNtutorial/getting_started.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/book/content/UNtutorial/getting_started.md b/docs/book/content/UNtutorial/getting_started.md index 33e5baa..5fc50d1 100644 --- a/docs/book/content/UNtutorial/getting_started.md +++ b/docs/book/content/UNtutorial/getting_started.md @@ -32,7 +32,7 @@ For the October 21-25, 2024 United Nations `OG-PHL` training in Manila, we will - Theory: ["Simple" 3-period-lived agent model](https://eapd-drb.github.io/OG-PHL/content/UNtutorial/3perOG.html) * - Tue. - Morning - - Review 3-period-lived-agent exercises (solutions in [this folder](https://github.com/EAPD-DRB/OG-PHL/blob/main/docs/book/content/UNtutorial/Solutions/3perOG/))
Review OG-Core and OG-PHL modules
Quick Git and GitHub workflow review + - Review 3-period-lived-agent exercises (solutions in [this folder](https://github.com/EAPD-DRB/OG-PHL/blob/main/docs/book/content/UNtutorial/solutions/3perOG/))
Review OG-Core and OG-PHL modules
Quick Git and GitHub workflow review * - - Afternoon - Running OG-PHL, inputs, outputs
Calibrating OG-PHL, current state, still to do From 480fb5734dd009e1b4904d7c04ba2c4d76b168ea Mon Sep 17 00:00:00 2001 From: Richard Evans Date: Wed, 23 Oct 2024 00:24:04 +0800 Subject: [PATCH 3/4] Black formatted solution files --- .../content/UNtutorial/solutions/3perOG/SS.py | 188 +++++++----- .../UNtutorial/solutions/3perOG/TPI.py | 276 +++++++++++------- .../UNtutorial/solutions/3perOG/execute.py | 163 +++++++---- 3 files changed, 388 insertions(+), 239 deletions(-) diff --git a/docs/book/content/UNtutorial/solutions/3perOG/SS.py b/docs/book/content/UNtutorial/solutions/3perOG/SS.py index 231bd87..7f4840d 100644 --- a/docs/book/content/UNtutorial/solutions/3perOG/SS.py +++ b/docs/book/content/UNtutorial/solutions/3perOG/SS.py @@ -1,4 +1,4 @@ -''' +""" ------------------------------------------------------------------------ This module contains the functions used to solve the steady state for the model with 3-period lived agents and exogenous labor from Chapter 5 @@ -18,7 +18,8 @@ get_b_errors() get_SS() ------------------------------------------------------------------------ -''' +""" + # Import packages import time import numpy as np @@ -32,7 +33,7 @@ def print_time(seconds, type): - ''' + """ -------------------------------------------------------------------- Takes a total amount of time in seconds and prints it in terms of more readable units (days, hours, minutes, seconds) @@ -53,34 +54,58 @@ def print_time(seconds, type): RETURNS: Nothing -------------------------------------------------------------------- - ''' + """ if seconds < 60: # seconds secs = round(seconds, 4) - print(type + ' computation time: ' + str(secs) + ' sec') + print(type + " computation time: " + str(secs) + " sec") elif seconds >= 60 and seconds < 3600: # minutes mins = int(seconds / 60) secs = round(((seconds / 60) - mins) * 60, 1) - print(type + ' computation time: ' + str(mins) + ' min, ' + - str(secs) + ' sec') + print( + type + + " computation time: " + + str(mins) + + " min, " + + str(secs) + + " sec" + ) elif seconds >= 3600 and seconds < 86400: # hours hrs = int(seconds / 3600) mins = int(((seconds / 3600) - hrs) * 60) secs = round(((seconds / 60) - hrs * 60 - mins) * 60, 1) - print(type + ' computation time: ' + str(hrs) + ' hrs, ' + - str(mins) + ' min, ' + str(secs) + ' sec') + print( + type + + " computation time: " + + str(hrs) + + " hrs, " + + str(mins) + + " min, " + + str(secs) + + " sec" + ) elif seconds >= 86400: # days days = int(seconds / 86400) hrs = int(((seconds / 86400) - days) * 24) mins = int(((seconds / 3600) - days * 24 - hrs) * 60) secs = round( - ((seconds / 60) - days * 24 * 60 - hrs * 60 - mins) * 60, 1) - print(type + ' computation time: ' + str(days) + ' days, ' + - str(hrs) + ' hrs, ' + str(mins) + ' min, ' + - str(secs) + ' sec') + ((seconds / 60) - days * 24 * 60 - hrs * 60 - mins) * 60, 1 + ) + print( + type + + " computation time: " + + str(days) + + " days, " + + str(hrs) + + " hrs, " + + str(mins) + + " min, " + + str(secs) + + " sec" + ) def get_cvec(r, w, bvec, nvec): - ''' + """ -------------------------------------------------------------------- Generates vector of lifetime steady-state consumptions given savings decisions, parameters, and the corresponding steady-state interest @@ -106,20 +131,22 @@ def get_cvec(r, w, bvec, nvec): RETURNS: cvec, c_cnstr -------------------------------------------------------------------- - ''' + """ b_s = np.append([0], bvec) b_sp1 = np.append(bvec, [0]) cvec = (1 + r) * b_s + w * nvec - b_sp1 if cvec.min() <= 0: - print('get_cvec() warning: distribution of savings and/or ' + - 'parameters created c<=0 for some agent(s)') + print( + "get_cvec() warning: distribution of savings and/or " + + "parameters created c<=0 for some agent(s)" + ) c_cnstr = cvec <= 0 return cvec, c_cnstr def get_L(nvec): - ''' + """ -------------------------------------------------------------------- Solve for aggregate labor L -------------------------------------------------------------------- @@ -135,14 +162,14 @@ def get_L(nvec): RETURNS: L -------------------------------------------------------------------- - ''' + """ L = nvec.sum() return L def get_K(barr): - ''' + """ -------------------------------------------------------------------- Solve for steady-state aggregate capital stock K or time path of aggregate capital stock K_t @@ -162,26 +189,30 @@ def get_K(barr): RETURNS: K, K_cnstr -------------------------------------------------------------------- - ''' + """ if barr.ndim == 1: # This is the steady-state case K = barr.sum() K_cnstr = K <= 0 if K_cnstr: - print('get_K() warning: distribution of savings and/or ' + - 'parameters created K<=0 for some agent(s)') + print( + "get_K() warning: distribution of savings and/or " + + "parameters created K<=0 for some agent(s)" + ) elif barr.ndim == 2: # This is the time path case K = barr.sum(axis=0) K_cnstr = K <= 0 if K.min() <= 0: - print('Aggregate capital constraint is violated K<=0 for ' + - 'some period in time path.') + print( + "Aggregate capital constraint is violated K<=0 for " + + "some period in time path." + ) return K, K_cnstr def get_w(params, K, L): - ''' + """ -------------------------------------------------------------------- Solve for steady-state wage w or time path of wages w_t -------------------------------------------------------------------- @@ -204,7 +235,7 @@ def get_w(params, K, L): RETURNS: w -------------------------------------------------------------------- - ''' + """ A, alpha = params w = (1 - alpha) * A * ((K / L) ** alpha) @@ -212,7 +243,7 @@ def get_w(params, K, L): def get_r(params, K, L): - ''' + """ -------------------------------------------------------------------- Solve for steady-state interest rate r or time path of interest rates r_t @@ -237,7 +268,7 @@ def get_r(params, K, L): RETURNS: r -------------------------------------------------------------------- - ''' + """ A, alpha, delta = params r = alpha * A * ((L / K) ** (1 - alpha)) - delta @@ -245,7 +276,7 @@ def get_r(params, K, L): def get_Y(params, K, L): - ''' + """ -------------------------------------------------------------------- Solve for steady-state aggregate output Y or time path of aggregate output Y_t @@ -270,15 +301,15 @@ def get_Y(params, K, L): RETURNS: Y -------------------------------------------------------------------- - ''' + """ A, alpha = params - Y = A * (K ** alpha) * (L ** (1 - alpha)) + Y = A * (K**alpha) * (L ** (1 - alpha)) return Y def get_C(carr): - ''' + """ -------------------------------------------------------------------- Solve for steady-state aggregate consumption C or time path of aggregate consumption C_t @@ -296,7 +327,7 @@ def get_C(carr): Returns: C -------------------------------------------------------------------- - ''' + """ if carr.ndim == 1: C = carr.sum() elif carr.ndim == 2: @@ -306,7 +337,7 @@ def get_C(carr): def feasible(params, bvec): - ''' + """ -------------------------------------------------------------------- Check whether a vector of steady-state savings is feasible in that it satisfies the nonnegativity constraints on consumption in every @@ -345,7 +376,7 @@ def feasible(params, bvec): RETURNS: b_cnstr, c_cnstr, K_cnstr -------------------------------------------------------------------- - ''' + """ nvec, A, alpha, delta = params L = get_L(nvec) K, K_cnstr = get_K(bvec) @@ -364,7 +395,7 @@ def feasible(params, bvec): def EulerSys(bvec, *args): - ''' + """ -------------------------------------------------------------------- Generates vector of all Euler errors that characterize optimal lifetime decisions @@ -408,11 +439,11 @@ def EulerSys(bvec, *args): RETURNS: b_errors -------------------------------------------------------------------- - ''' + """ beta, sigma, nvec, L, A, alpha, delta, EulDiff = args K, K_cnstr = get_K(bvec) if K_cnstr: - b_err_vec = 1000. * np.ones(nvec.shape[0] - 1) + b_err_vec = 1000.0 * np.ones(nvec.shape[0] - 1) else: r_params = (A, alpha, delta) r = get_r(r_params, K, L) @@ -420,14 +451,13 @@ def EulerSys(bvec, *args): w = get_w(w_params, K, L) cvec, c_cnstr = get_cvec(r, w, bvec, nvec) b_err_params = (beta, sigma) - b_err_vec = get_b_errors(b_err_params, r, cvec, c_cnstr, - EulDiff) + b_err_vec = get_b_errors(b_err_params, r, cvec, c_cnstr, EulDiff) return b_err_vec def get_b_errors(params, r, cvec, c_cnstr, diff): - ''' + """ -------------------------------------------------------------------- Generates vector of dynamic Euler errors that characterize the optimal lifetime savings decision. Because this function is used for @@ -459,26 +489,26 @@ def get_b_errors(params, r, cvec, c_cnstr, diff): RETURNS: b_errors -------------------------------------------------------------------- - ''' + """ beta, sigma = params # Make each negative consumption artifically positive - cvec[c_cnstr] = 9999. + cvec[c_cnstr] = 9999.0 mu_c = cvec[:-1] ** (-sigma) mu_cp1 = cvec[1:] ** (-sigma) if diff: b_errors = (beta * (1 + r) * mu_cp1) - mu_c - b_errors[c_cnstr[:-1]] = 9999. - b_errors[c_cnstr[1:]] = 9999. + b_errors[c_cnstr[:-1]] = 9999.0 + b_errors[c_cnstr[1:]] = 9999.0 else: b_errors = ((beta * (1 + r) * mu_cp1) / mu_c) - 1 - b_errors[c_cnstr[:-1]] = 9999. / 100 - b_errors[c_cnstr[1:]] = 9999. / 100 + b_errors[c_cnstr[:-1]] = 9999.0 / 100 + b_errors[c_cnstr[1:]] = 9999.0 / 100 return b_errors def get_SS(params, bvec_guess, graphs): - ''' + """ -------------------------------------------------------------------- Solve for the steady-state solution of the 3-period-lived agent OG model with exogenous labor supply @@ -539,21 +569,24 @@ def get_SS(params, bvec_guess, graphs): RETURNS: ss_output -------------------------------------------------------------------- - ''' + """ start_time = time.time() beta, sigma, nvec, L, A, alpha, delta, SS_tol, EulDiff = params f_params = (nvec, A, alpha, delta) b1_cnstr, c1_cnstr, K1_cnstr = feasible(f_params, bvec_guess) if K1_cnstr is True or c1_cnstr.max() is True: - err = ("Initial guess problem: " + - "One or more constraints not satisfied.") + err = ( + "Initial guess problem: " + + "One or more constraints not satisfied." + ) print("K1_cnstr: ", K1_cnstr) print("c1_cnstr: ", c1_cnstr) raise RuntimeError(err) else: eul_args = (beta, sigma, nvec, L, A, alpha, delta, EulDiff) - results_bss = opt.root(EulerSys, bvec_guess, args=(eul_args), - tol=SS_tol) + results_bss = opt.root( + EulerSys, bvec_guess, args=(eul_args), tol=SS_tol + ) b_ss = results_bss.x # Generate other steady-state values and Euler equations @@ -567,27 +600,33 @@ def get_SS(params, bvec_guess, graphs): Y_ss = get_Y(Y_params, K_ss, L) C_ss = get_C(c_ss) b_err_params = (beta, sigma) - EulErr_ss = get_b_errors( - b_err_params, r_ss, c_ss, c_cnstr, EulDiff) + EulErr_ss = get_b_errors(b_err_params, r_ss, c_ss, c_cnstr, EulDiff) RCerr_ss = Y_ss - C_ss - delta * K_ss ss_time = time.time() - start_time ss_output = { - 'b_ss': b_ss, 'c_ss': c_ss, 'w_ss': w_ss, 'r_ss': r_ss, - 'K_ss': K_ss, 'Y_ss': Y_ss, 'C_ss': C_ss, - 'EulErr_ss': EulErr_ss, 'RCerr_ss': RCerr_ss, - 'ss_time': ss_time} - print('b_ss is: ', b_ss) - print('c_ss is: ', c_ss) - print('Euler errors are: ', EulErr_ss) - print('Resource constraint error is: ', RCerr_ss) + "b_ss": b_ss, + "c_ss": c_ss, + "w_ss": w_ss, + "r_ss": r_ss, + "K_ss": K_ss, + "Y_ss": Y_ss, + "C_ss": C_ss, + "EulErr_ss": EulErr_ss, + "RCerr_ss": RCerr_ss, + "ss_time": ss_time, + } + print("b_ss is: ", b_ss) + print("c_ss is: ", c_ss) + print("Euler errors are: ", EulErr_ss) + print("Resource constraint error is: ", RCerr_ss) # Print SS computation time - print_time(ss_time, 'SS') + print_time(ss_time, "SS") if graphs: - ''' + """ ---------------------------------------------------------------- cur_path = string, path name of current directory output_fldr = string, folder in current path to save files @@ -596,7 +635,7 @@ def get_SS(params, bvec_guess, graphs): S = integer >= 3, number of periods in a life age_pers = (S,) vector, ages from 1 to S ---------------------------------------------------------------- - ''' + """ # Create directory if images directory does not already exist cur_path = os.path.split(os.path.abspath(__file__))[0] output_fldr = "images" @@ -608,19 +647,18 @@ def get_SS(params, bvec_guess, graphs): S = nvec.shape[0] age_pers = np.arange(1, S + 1) fig, ax = plt.subplots() - plt.plot(age_pers, c_ss, marker='D', label='Consumption') - plt.plot(age_pers, np.hstack((0, b_ss)), marker='D', - label='Savings') + plt.plot(age_pers, c_ss, marker="D", label="Consumption") + plt.plot(age_pers, np.hstack((0, b_ss)), marker="D", label="Savings") # for the minor ticks, use no labels; default NullFormatter minorLocator = MultipleLocator(1) ax.xaxis.set_minor_locator(minorLocator) - plt.grid(visible=True, which='major', color='0.65', linestyle='-') - plt.title('Steady-state consumption and savings', fontsize=20) - plt.xlabel(r'Age $s$') - plt.ylabel(r'Units of consumption') + plt.grid(visible=True, which="major", color="0.65", linestyle="-") + plt.title("Steady-state consumption and savings", fontsize=20) + plt.xlabel(r"Age $s$") + plt.ylabel(r"Units of consumption") plt.xlim((0.8, S + 0.2)) plt.ylim((-0.02, 1.15 * (c_ss.max()))) - plt.legend(loc='center right') + plt.legend(loc="center right") output_path = os.path.join(output_dir, "SS_bc") plt.savefig(output_path) # plt.show() diff --git a/docs/book/content/UNtutorial/solutions/3perOG/TPI.py b/docs/book/content/UNtutorial/solutions/3perOG/TPI.py index 7f53bc3..038336a 100644 --- a/docs/book/content/UNtutorial/solutions/3perOG/TPI.py +++ b/docs/book/content/UNtutorial/solutions/3perOG/TPI.py @@ -1,4 +1,4 @@ -''' +""" ------------------------------------------------------------------------ This module contains the functions used to solve the time path iteration non-steady-state equilibrium for the model with 3-period lived agents @@ -25,7 +25,8 @@ get_b_errors() get_SS() ------------------------------------------------------------------------ -''' +""" + # Import Packages import time import numpy as np @@ -38,15 +39,15 @@ import sys import os -''' +""" ------------------------------------------------------------------------ Functions ------------------------------------------------------------------------ -''' +""" def get_path(x1, xT, T, spec): - ''' + """ -------------------------------------------------------------------- This function generates a path from point x1 to point xT such that that the path x is a linear or quadratic function of time t. @@ -78,7 +79,7 @@ def get_path(x1, xT, T, spec): RETURNS: xpath -------------------------------------------------------------------- - ''' + """ if spec == "linear": xpath = np.linspace(x1, xT, T) elif spec == "quadratic": @@ -91,7 +92,7 @@ def get_path(x1, xT, T, spec): def get_cvec_lf(rpath, wpath, nvec, bvec): - ''' + """ -------------------------------------------------------------------- Generates vector of remaining lifetime consumptions from individual savings, and the time path of interest rates and the real wages, @@ -117,19 +118,21 @@ def get_cvec_lf(rpath, wpath, nvec, bvec): RETURNS: cvec, c_cnstr -------------------------------------------------------------------- - ''' + """ b_s = bvec b_sp1 = np.append(bvec[1:], [0]) cvec = (1 + rpath) * b_s + wpath * nvec - b_sp1 if cvec.min() <= 0: - print('get_cvec_lf() warning: distribution of savings and/or ' + - 'parameters created c<=0 for some agent(s)') + print( + "get_cvec_lf() warning: distribution of savings and/or " + + "parameters created c<=0 for some agent(s)" + ) c_cnstr = cvec <= 0 return cvec, c_cnstr def LfEulerSys(bvec, *args): - ''' + """ -------------------------------------------------------------------- Generates vector of all Euler errors for a given bvec, which errors characterize all optimal lifetime decisions, where p is an integer @@ -165,19 +168,19 @@ def LfEulerSys(bvec, *args): RETURNS: b_err_vec -------------------------------------------------------------------- - ''' + """ beta, sigma, beg_wealth, nvec, rpath, wpath, EulDiff = args bvec2 = np.append(beg_wealth, bvec) cvec, c_cnstr = get_cvec_lf(rpath, wpath, nvec, bvec2) b_err_params = (beta, sigma) - b_err_vec = ss.get_b_errors(b_err_params, rpath[1:], cvec, - c_cnstr, EulDiff) + b_err_vec = ss.get_b_errors( + b_err_params, rpath[1:], cvec, c_cnstr, EulDiff + ) return b_err_vec -def paths_life(params, beg_age, beg_wealth, nvec, rpath, wpath, - b_init): - ''' +def paths_life(params, beg_age, beg_wealth, nvec, rpath, wpath, b_init): + """ -------------------------------------------------------------------- Solve for the remaining lifetime savings decisions of an individual who enters the model at age beg_age, with corresponding initial @@ -227,7 +230,7 @@ def paths_life(params, beg_age, beg_wealth, nvec, rpath, wpath, RETURNS: bpath, cpath, b_err_vec -------------------------------------------------------------------- - ''' + """ S, beta, sigma, TPI_tol, EulDiff = params p = int(S - beg_age + 1) if beg_age == 1 and beg_wealth != 0: @@ -240,18 +243,19 @@ def paths_life(params, beg_age, beg_wealth, nvec, rpath, wpath, sys.exit("Beginning age and length of nvec do not match.") b_guess = 1.01 * b_init eullf_objs = (beta, sigma, beg_wealth, nvec, rpath, wpath, EulDiff) - bpath = opt.fsolve(LfEulerSys, b_guess, args=(eullf_objs), - xtol=TPI_tol) - cpath, c_cnstr = get_cvec_lf(rpath, wpath, nvec, - np.append(beg_wealth, bpath)) + bpath = opt.fsolve(LfEulerSys, b_guess, args=(eullf_objs), xtol=TPI_tol) + cpath, c_cnstr = get_cvec_lf( + rpath, wpath, nvec, np.append(beg_wealth, bpath) + ) b_err_params = (beta, sigma) - b_err_vec = ss.get_b_errors(b_err_params, rpath[1:], cpath, - c_cnstr, EulDiff) + b_err_vec = ss.get_b_errors( + b_err_params, rpath[1:], cpath, c_cnstr, EulDiff + ) return bpath, cpath, b_err_vec def get_cbepath(params, rpath, wpath): - ''' + """ -------------------------------------------------------------------- Generates matrices for the time path of the distribution of individual savings, individual consumption, and the Euler errors @@ -302,54 +306,68 @@ def get_cbepath(params, rpath, wpath): RETURNS: cpath, bpath, EulErrPath -------------------------------------------------------------------- - ''' + """ S, T, beta, sigma, nvec, bvec1, b_ss, TPI_tol, EulDiff = params cpath = np.zeros((S, T + S - 2)) - bpath = np.append(bvec1.reshape((S - 1, 1)), - np.zeros((S - 1, T + S - 3)), axis=1) + bpath = np.append( + bvec1.reshape((S - 1, 1)), np.zeros((S - 1, T + S - 3)), axis=1 + ) EulErrPath = np.zeros((S - 1, T + S - 2)) # Solve the incomplete remaining lifetime decisions of agents alive # in period t=1 but not born in period t=1 - cpath[S - 1, 0] = ((1 + rpath[0]) * bvec1[S - 2] + - wpath[0] * nvec[S - 1]) + cpath[S - 1, 0] = (1 + rpath[0]) * bvec1[S - 2] + wpath[0] * nvec[S - 1] pl_params = (S, beta, sigma, TPI_tol, EulDiff) for p in range(2, S): - b_guess = np.diagonal(bpath[S - p:, :p - 1]) + b_guess = np.diagonal(bpath[S - p :, : p - 1]) bveclf, cveclf, b_err_veclf = paths_life( - pl_params, S - p + 1, bvec1[S - p - 1], nvec[-p:], - rpath[:p], wpath[:p], b_guess) + pl_params, + S - p + 1, + bvec1[S - p - 1], + nvec[-p:], + rpath[:p], + wpath[:p], + b_guess, + ) # Insert the vector lifetime solutions diagonally (twist donut) # into the cpath, bpath, and EulErrPath matrices DiagMaskb = np.eye(p - 1, dtype=bool) DiagMaskc = np.eye(p, dtype=bool) - bpath[S - p:, 1:p] = DiagMaskb * bveclf + bpath[S - p:, 1:p] - cpath[S - p:, :p] = DiagMaskc * cveclf + cpath[S - p:, :p] - EulErrPath[S - p:, 1:p] = (DiagMaskb * b_err_veclf + - EulErrPath[S - p:, 1:p]) + bpath[S - p :, 1:p] = DiagMaskb * bveclf + bpath[S - p :, 1:p] + cpath[S - p :, :p] = DiagMaskc * cveclf + cpath[S - p :, :p] + EulErrPath[S - p :, 1:p] = ( + DiagMaskb * b_err_veclf + EulErrPath[S - p :, 1:p] + ) # Solve for complete lifetime decisions of agents born in periods # 1 to T and insert the vector lifetime solutions diagonally (twist # donut) into the cpath, bpath, and EulErrPath matrices DiagMaskb = np.eye(S - 1, dtype=bool) DiagMaskc = np.eye(S, dtype=bool) for t in range(1, T): # Go from periods 1 to T-1 - b_guess = np.diagonal(bpath[:, t - 1:t + S - 2]) + b_guess = np.diagonal(bpath[:, t - 1 : t + S - 2]) bveclf, cveclf, b_err_veclf = paths_life( - pl_params, 1, 0, nvec, rpath[t - 1:t + S - 1], - wpath[t - 1:t + S - 1], b_guess) + pl_params, + 1, + 0, + nvec, + rpath[t - 1 : t + S - 1], + wpath[t - 1 : t + S - 1], + b_guess, + ) # Insert the vector lifetime solutions diagonally (twist donut) # into the cpath, bpath, and EulErrPath matrices - bpath[:, t:t + S - 1] = (DiagMaskb * bveclf + - bpath[:, t:t + S - 1]) - cpath[:, t - 1:t + S - 1] = (DiagMaskc * cveclf + - cpath[:, t - 1:t + S - 1]) - EulErrPath[:, t:t + S - 1] = (DiagMaskb * b_err_veclf + - EulErrPath[:, t:t + S - 1]) + bpath[:, t : t + S - 1] = DiagMaskb * bveclf + bpath[:, t : t + S - 1] + cpath[:, t - 1 : t + S - 1] = ( + DiagMaskc * cveclf + cpath[:, t - 1 : t + S - 1] + ) + EulErrPath[:, t : t + S - 1] = ( + DiagMaskb * b_err_veclf + EulErrPath[:, t : t + S - 1] + ) return cpath, bpath, EulErrPath def get_TPI(params, bvec1, graphs): - ''' + """ -------------------------------------------------------------------- Generates steady-state time path for all endogenous objects from initial state (K1, Gamma1) to the steady state. @@ -445,10 +463,27 @@ def get_TPI(params, bvec1, graphs): RETURNS: tpi_output -------------------------------------------------------------------- - ''' + """ start_time = time.time() - (S, T, beta, sigma, nvec, L, A, alpha, delta, b_ss, K_ss, C_ss, - maxiter_TPI, mindist_TPI, xi, TPI_tol, EulDiff) = params + ( + S, + T, + beta, + sigma, + nvec, + L, + A, + alpha, + delta, + b_ss, + K_ss, + C_ss, + maxiter_TPI, + mindist_TPI, + xi, + TPI_tol, + EulDiff, + ) = params K1, K1_cnstr = ss.get_K(bvec1) # Create time paths for K and L @@ -458,12 +493,11 @@ def get_TPI(params, bvec1, graphs): Lpath = L * np.ones(T + S - 2) iter_TPI = int(0) - dist_TPI = 10. + dist_TPI = 10.0 Kpath_new = Kpath_init.copy() r_params = (A, alpha, delta) w_params = (A, alpha) - cbe_params = (S, T, beta, sigma, nvec, bvec1, b_ss, TPI_tol, - EulDiff) + cbe_params = (S, T, beta, sigma, nvec, bvec1, b_ss, TPI_tol, EulDiff) while (iter_TPI < maxiter_TPI) and (dist_TPI >= mindist_TPI): iter_TPI += 1 @@ -474,42 +508,55 @@ def get_TPI(params, bvec1, graphs): Kpath_new = np.zeros(T + S - 2) Kpath_new[:T], Kpath_cnstr = ss.get_K(bpath[:, :T]) Kpath_new[T:] = K_ss * np.ones(S - 2) - Kpath_cnstr = np.append(Kpath_cnstr, - np.zeros(S - 2, dtype=bool)) + Kpath_cnstr = np.append(Kpath_cnstr, np.zeros(S - 2, dtype=bool)) Kpath_new[Kpath_cnstr] = 0.1 # Check the distance of Kpath_new1 dist_TPI = ((Kpath_new[1:T] - Kpath_init[1:T]) ** 2).sum() # dist_TPI = np.absolute((Kpath_new[1:T] - Kpath_init[1:T]) / # Kpath_init[1:T]).max() - print('iter: ', iter_TPI, ', dist: ', dist_TPI, - ',max Eul err: ', np.absolute(EulErrPath).max()) + print( + "iter: ", + iter_TPI, + ", dist: ", + dist_TPI, + ",max Eul err: ", + np.absolute(EulErrPath).max(), + ) if iter_TPI == maxiter_TPI and dist_TPI > mindist_TPI: - print('TPI reached maxiter and did not converge.') + print("TPI reached maxiter and did not converge.") elif iter_TPI == maxiter_TPI and dist_TPI <= mindist_TPI: - print('TPI converged in the last iteration. ' + - 'Should probably increase maxiter_TPI.') + print( + "TPI converged in the last iteration. " + + "Should probably increase maxiter_TPI." + ) Kpath = Kpath_new Y_params = (A, alpha) Ypath = ss.get_Y(Y_params, Kpath, Lpath) Cpath = np.zeros(T + S - 2) - Cpath[:T - 1] = ss.get_C(cpath[:, :T - 1]) - Cpath[T - 1:] = C_ss * np.ones(S - 1) - RCerrPath = (Ypath[:-1] - Cpath[:-1] - Kpath[1:] + - (1 - delta) * Kpath[:-1]) + Cpath[: T - 1] = ss.get_C(cpath[:, : T - 1]) + Cpath[T - 1 :] = C_ss * np.ones(S - 1) + RCerrPath = Ypath[:-1] - Cpath[:-1] - Kpath[1:] + (1 - delta) * Kpath[:-1] tpi_time = time.time() - start_time tpi_output = { - 'bpath': bpath, 'cpath': cpath, 'wpath': wpath, 'rpath': rpath, - 'Kpath': Kpath, 'Ypath': Ypath, 'Cpath': Cpath, - 'EulErrPath': EulErrPath, 'RCerrPath': RCerrPath, - 'tpi_time': tpi_time} + "bpath": bpath, + "cpath": cpath, + "wpath": wpath, + "rpath": rpath, + "Kpath": Kpath, + "Ypath": Ypath, + "Cpath": Cpath, + "EulErrPath": EulErrPath, + "RCerrPath": RCerrPath, + "tpi_time": tpi_time, + } # Print TPI computation time - ss.print_time(tpi_time, 'TPI') + ss.print_time(tpi_time, "TPI") if graphs: - ''' + """ ---------------------------------------------------------------- cur_path = string, path name of current directory output_fldr = string, folder in current path to save files @@ -529,7 +576,7 @@ def get_TPI(params, bvec1, graphs): smatc = (3, 17) matrix, ages for all consumption decisions ages (S) and time periods (T-1) ---------------------------------------------------------------- - ''' + """ # Create directory if images directory does not already exist cur_path = os.path.split(os.path.abspath(__file__))[0] output_fldr = "images" @@ -541,65 +588,65 @@ def get_TPI(params, bvec1, graphs): tvec = np.linspace(1, T + S - 2, T + S - 2) minorLocator = MultipleLocator(1) fig, ax = plt.subplots() - plt.plot(tvec, Kpath, marker='D') + plt.plot(tvec, Kpath, marker="D") # for the minor ticks, use no labels; default NullFormatter ax.xaxis.set_minor_locator(minorLocator) - plt.grid(visible=True, which='major', color='0.65', linestyle='-') - plt.title('Time path for aggregate capital stock K') - plt.xlabel(r'Period $t$') - plt.ylabel(r'Aggregate capital $K_{t}$') + plt.grid(visible=True, which="major", color="0.65", linestyle="-") + plt.title("Time path for aggregate capital stock K") + plt.xlabel(r"Period $t$") + plt.ylabel(r"Aggregate capital $K_{t}$") output_path = os.path.join(output_dir, "Kpath") plt.savefig(output_path) # plt.show() # Plot time path of aggregate output (GDP) fig, ax = plt.subplots() - plt.plot(tvec, Ypath, marker='D') + plt.plot(tvec, Ypath, marker="D") # for the minor ticks, use no labels; default NullFormatter ax.xaxis.set_minor_locator(minorLocator) - plt.grid(visible=True, which='major', color='0.65', linestyle='-') - plt.title('Time path for aggregate output (GDP) Y') - plt.xlabel(r'Period $t$') - plt.ylabel(r'Aggregate output $Y_{t}$') + plt.grid(visible=True, which="major", color="0.65", linestyle="-") + plt.title("Time path for aggregate output (GDP) Y") + plt.xlabel(r"Period $t$") + plt.ylabel(r"Aggregate output $Y_{t}$") output_path = os.path.join(output_dir, "Ypath") plt.savefig(output_path) # plt.show() # Plot time path of aggregate consumption fig, ax = plt.subplots() - plt.plot(tvec, Cpath, marker='D') + plt.plot(tvec, Cpath, marker="D") # for the minor ticks, use no labels; default NullFormatter ax.xaxis.set_minor_locator(minorLocator) - plt.grid(visible=True, which='major', color='0.65', linestyle='-') - plt.title('Time path for aggregate consumption C') - plt.xlabel(r'Period $t$') - plt.ylabel(r'Aggregate consumption $C_{t}$') + plt.grid(visible=True, which="major", color="0.65", linestyle="-") + plt.title("Time path for aggregate consumption C") + plt.xlabel(r"Period $t$") + plt.ylabel(r"Aggregate consumption $C_{t}$") output_path = os.path.join(output_dir, "C_aggr_path") plt.savefig(output_path) # plt.show() # Plot time path of real wage fig, ax = plt.subplots() - plt.plot(tvec, wpath, marker='D') + plt.plot(tvec, wpath, marker="D") # for the minor ticks, use no labels; default NullFormatter ax.xaxis.set_minor_locator(minorLocator) - plt.grid(visible=True, which='major', color='0.65', linestyle='-') - plt.title('Time path for real wage w') - plt.xlabel(r'Period $t$') - plt.ylabel(r'Real wage $w_{t}$') + plt.grid(visible=True, which="major", color="0.65", linestyle="-") + plt.title("Time path for real wage w") + plt.xlabel(r"Period $t$") + plt.ylabel(r"Real wage $w_{t}$") output_path = os.path.join(output_dir, "wpath") plt.savefig(output_path) # plt.show() # Plot time path of real interest rate fig, ax = plt.subplots() - plt.plot(tvec, rpath, marker='D') + plt.plot(tvec, rpath, marker="D") # for the minor ticks, use no labels; default NullFormatter ax.xaxis.set_minor_locator(minorLocator) - plt.grid(visible=True, which='major', color='0.65', linestyle='-') - plt.title('Time path for real interest rate r') - plt.xlabel(r'Period $t$') - plt.ylabel(r'Real interest rate $r_{t}$') + plt.grid(visible=True, which="major", color="0.65", linestyle="-") + plt.title("Time path for real interest rate r") + plt.xlabel(r"Period $t$") + plt.ylabel(r"Real interest rate $r_{t}$") output_path = os.path.join(output_dir, "rpath") plt.savefig(output_path) # plt.show() @@ -608,14 +655,20 @@ def get_TPI(params, bvec1, graphs): tgridT = np.linspace(1, T, T) sgrid2 = np.linspace(2, S, S - 1) tmatb, smatb = np.meshgrid(tgridT, sgrid2) - cmap_bp = matplotlib.cm.get_cmap('summer') + cmap_bp = matplotlib.cm.get_cmap("summer") fig, ax = plt.subplots(subplot_kw={"projection": "3d"}) - ax.set_xlabel(r'period-$t$') - ax.set_ylabel(r'age-$s$') - ax.set_zlabel(r'individual savings $b_{s,t}$') + ax.set_xlabel(r"period-$t$") + ax.set_ylabel(r"age-$s$") + ax.set_zlabel(r"individual savings $b_{s,t}$") strideval = max(int(1), int(round(S / 10))) - ax.plot_surface(tmatb, smatb, bpath[:, :T], rstride=strideval, - cstride=strideval, cmap=cmap_bp) + ax.plot_surface( + tmatb, + smatb, + bpath[:, :T], + rstride=strideval, + cstride=strideval, + cmap=cmap_bp, + ) output_path = os.path.join(output_dir, "bpath") plt.savefig(output_path) # plt.show() @@ -624,15 +677,20 @@ def get_TPI(params, bvec1, graphs): tgridTm1 = np.linspace(1, T - 1, T - 1) sgrid = np.linspace(1, S, S) tmatc, smatc = np.meshgrid(tgridTm1, sgrid) - cmap_cp = matplotlib.cm.get_cmap('summer') + cmap_cp = matplotlib.cm.get_cmap("summer") fig, ax = plt.subplots(subplot_kw={"projection": "3d"}) - ax.set_xlabel(r'period-$t$') - ax.set_ylabel(r'age-$s$') - ax.set_zlabel(r'individual consumption $c_{s,t}$') + ax.set_xlabel(r"period-$t$") + ax.set_ylabel(r"age-$s$") + ax.set_zlabel(r"individual consumption $c_{s,t}$") strideval = max(int(1), int(round(S / 10))) - ax.plot_surface(tmatc, smatc, cpath[:, :T - 1], - rstride=strideval, cstride=strideval, - cmap=cmap_cp) + ax.plot_surface( + tmatc, + smatc, + cpath[:, : T - 1], + rstride=strideval, + cstride=strideval, + cmap=cmap_cp, + ) output_path = os.path.join(output_dir, "cpath") plt.savefig(output_path) # plt.show() diff --git a/docs/book/content/UNtutorial/solutions/3perOG/execute.py b/docs/book/content/UNtutorial/solutions/3perOG/execute.py index 41afb43..6efa79d 100755 --- a/docs/book/content/UNtutorial/solutions/3perOG/execute.py +++ b/docs/book/content/UNtutorial/solutions/3perOG/execute.py @@ -1,4 +1,4 @@ -''' +""" ------------------------------------------------------------------------ This program runs the steady state solver as well as the time path iteration solution for the model with 3-period lived agents and @@ -14,13 +14,14 @@ ss.get_K tpi.get_TPI() ------------------------------------------------------------------------ -''' +""" + # Import packages import numpy as np import SS as ss import TPI as tpi -''' +""" ------------------------------------------------------------------------ Declare parameters ------------------------------------------------------------------------ @@ -49,18 +50,18 @@ beta*(1+r)*u'(c2) - u'(c1), =False if want ratio version [beta*(1+r)*u'(c2)]/[u'(c1)] - 1 ------------------------------------------------------------------------ -''' +""" # Household parameters S = int(3) -beta_annual = .96 -beta = beta_annual ** 20 +beta_annual = 0.96 +beta = beta_annual**20 sigma = 3.0 nvec = np.array([1.0, 1.0, 0.2]) L = nvec.sum() # Firm parameters A = 1.0 -alpha = .35 -delta_annual = .05 +alpha = 0.35 +delta_annual = 0.05 delta = 1 - ((1 - delta_annual) ** 20) # SS parameters SS_tol = 1e-13 @@ -76,7 +77,7 @@ # Overall parameters EulDiff = False -''' +""" ------------------------------------------------------------------------ Check feasibility ------------------------------------------------------------------------ @@ -92,28 +93,28 @@ bvec_guess3 = (2,) vector, guess for steady-state bvec (b1, b2) ------------------------------------------------------------------------ -''' +""" f_params = (nvec, A, alpha, delta) bvec_guess1 = np.array([1.0, 1.2]) b_cnstr, c_cnstr, K_cnstr = ss.feasible(f_params, bvec_guess1) -print('bvec_guess1', bvec_guess1) -print('c_cnstr', c_cnstr) -print('K_cnstr', K_cnstr) +print("bvec_guess1", bvec_guess1) +print("c_cnstr", c_cnstr) +print("K_cnstr", K_cnstr) bvec_guess2 = np.array([0.06, -0.001]) b_cnstr, c_cnstr, K_cnstr = ss.feasible(f_params, bvec_guess2) -print('bvec_guess2', bvec_guess2) -print('c_cnstr', c_cnstr) -print('K_cnstr', K_cnstr) +print("bvec_guess2", bvec_guess2) +print("c_cnstr", c_cnstr) +print("K_cnstr", K_cnstr) bvec_guess3 = np.array([0.1, 0.1]) b_cnstr, c_cnstr, K_cnstr = ss.feasible(f_params, bvec_guess3) -print('bvec_guess3', bvec_guess3) -print('c_cnstr', c_cnstr) -print('K_cnstr', K_cnstr) +print("bvec_guess3", bvec_guess3) +print("c_cnstr", c_cnstr) +print("K_cnstr", K_cnstr) -''' +""" ------------------------------------------------------------------------ Run the steady-state solution ------------------------------------------------------------------------ @@ -130,42 +131,76 @@ ss_output = length 10 dictionary, {b_ss, c_ss, w_ss, r_ss, K_ss, Y_ss, C_ss, EulErr_ss, RCerr_ss, ss_time} ------------------------------------------------------------------------ -''' -print('BEGIN EQUILIBRIUM STEADY-STATE COMPUTATION') +""" +print("BEGIN EQUILIBRIUM STEADY-STATE COMPUTATION") bvec_guess = np.array([0.1, 0.1]) f_params = (nvec, A, alpha, delta) b_cnstr, c_cnstr, K_cnstr = ss.feasible(f_params, bvec_guess) if not K_cnstr and not c_cnstr.max(): - ss_params = (beta, sigma, nvec, L, A, alpha, delta, SS_tol, - EulDiff) + ss_params = (beta, sigma, nvec, L, A, alpha, delta, SS_tol, EulDiff) ss_output = ss.get_SS(ss_params, bvec_guess, SS_graphs) beta2 = 0.55 - ss_params2 = (beta2, sigma, nvec, L, A, alpha, delta, SS_tol, - EulDiff) + ss_params2 = (beta2, sigma, nvec, L, A, alpha, delta, SS_tol, EulDiff) ss_output2 = ss.get_SS(ss_params2, bvec_guess, False) - print('c2_ss = ', ss_output2['c_ss'], ', diff is ', - ss_output2['c_ss'] - ss_output['c_ss']) - print('b2_ss = ', ss_output2['b_ss'], ', diff is ', - ss_output2['b_ss'] - ss_output['b_ss']) - print('w2_ss = ', ss_output2['w_ss'], ', diff is ', - ss_output2['w_ss'] - ss_output['w_ss']) - print('r2_ss = ', ss_output2['r_ss'], ', diff is ', - ss_output2['r_ss'] - ss_output['r_ss']) - print('K2_ss = ', ss_output2['K_ss'], ', diff is ', - ss_output2['K_ss'] - ss_output['K_ss']) - print('Y2_ss = ', ss_output2['Y_ss'], ', diff is ', - ss_output2['Y_ss'] - ss_output['Y_ss']) - print('C2_ss = ', ss_output2['C_ss'], ', diff is ', - ss_output2['C_ss'] - ss_output['C_ss']) - print('EulErr2_ss = ', ss_output2['EulErr_ss'], ', diff is ', - ss_output2['EulErr_ss'] - ss_output['EulErr_ss']) - print('RCerr2_ss = ', ss_output2['RCerr_ss'], ', diff is ', - ss_output2['RCerr_ss'] - ss_output['RCerr_ss']) + print( + "c2_ss = ", + ss_output2["c_ss"], + ", diff is ", + ss_output2["c_ss"] - ss_output["c_ss"], + ) + print( + "b2_ss = ", + ss_output2["b_ss"], + ", diff is ", + ss_output2["b_ss"] - ss_output["b_ss"], + ) + print( + "w2_ss = ", + ss_output2["w_ss"], + ", diff is ", + ss_output2["w_ss"] - ss_output["w_ss"], + ) + print( + "r2_ss = ", + ss_output2["r_ss"], + ", diff is ", + ss_output2["r_ss"] - ss_output["r_ss"], + ) + print( + "K2_ss = ", + ss_output2["K_ss"], + ", diff is ", + ss_output2["K_ss"] - ss_output["K_ss"], + ) + print( + "Y2_ss = ", + ss_output2["Y_ss"], + ", diff is ", + ss_output2["Y_ss"] - ss_output["Y_ss"], + ) + print( + "C2_ss = ", + ss_output2["C_ss"], + ", diff is ", + ss_output2["C_ss"] - ss_output["C_ss"], + ) + print( + "EulErr2_ss = ", + ss_output2["EulErr_ss"], + ", diff is ", + ss_output2["EulErr_ss"] - ss_output["EulErr_ss"], + ) + print( + "RCerr2_ss = ", + ss_output2["RCerr_ss"], + ", diff is ", + ss_output2["RCerr_ss"] - ss_output["RCerr_ss"], + ) else: print("Initial guess for SS bvec does not satisfy K>0 or c_s>0.") -''' +""" ------------------------------------------------------------------------ Run the time path iteration (TPI) solution ------------------------------------------------------------------------ @@ -181,20 +216,38 @@ tpi_output = length 10 dictionary, {bpath, cpath, wpath, rpath, Kpath, Ypath, Cpath, EulErrPath, RCerrPath, tpi_time} ------------------------------------------------------------------------ -''' +""" if TPI_solve: - print('BEGIN EQUILIBRIUM TIME PATH COMPUTATION') - b_ss = ss_output['b_ss'] - K_ss = ss_output['K_ss'] - C_ss = ss_output['C_ss'] + print("BEGIN EQUILIBRIUM TIME PATH COMPUTATION") + b_ss = ss_output["b_ss"] + K_ss = ss_output["K_ss"] + C_ss = ss_output["C_ss"] bvec1 = np.array([0.8 * b_ss[0], 1.1 * b_ss[1]]) # Make sure init. period distribution is feasible in terms of K K1, K_constr_tpi1 = ss.get_K(bvec1) if K_constr_tpi1: - print('Initial savings distribution is not feasible because ' + - 'K1<=0. Some element(s) of bvec1 must increase.') + print( + "Initial savings distribution is not feasible because " + + "K1<=0. Some element(s) of bvec1 must increase." + ) else: - tpi_params = (S, T, beta, sigma, nvec, L, A, alpha, delta, b_ss, - K_ss, C_ss, maxiter_TPI, mindist_TPI, xi, TPI_tol, - EulDiff) + tpi_params = ( + S, + T, + beta, + sigma, + nvec, + L, + A, + alpha, + delta, + b_ss, + K_ss, + C_ss, + maxiter_TPI, + mindist_TPI, + xi, + TPI_tol, + EulDiff, + ) tpi_output = tpi.get_TPI(tpi_params, bvec1, TPI_graphs) From c275648b92f515285ccc3667acd76998a3f11ad6 Mon Sep 17 00:00:00 2001 From: Richard Evans Date: Wed, 23 Oct 2024 00:32:56 +0800 Subject: [PATCH 4/4] Added ipykernel to environment.yml --- environment.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/environment.yml b/environment.yml index ccade2e..1ef04a2 100644 --- a/environment.yml +++ b/environment.yml @@ -29,6 +29,7 @@ dependencies: - linearmodels - black - jupyter +- ipykernel - ipython - pip - pip: