#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sun Mar 21 08:40:21 2021
@author: A.Goumilevski
"""
import os
import numpy as np
import pandas as pd
from time import time
import scipy.sparse as sparse
from scipy.optimize import fmin
from scipy.linalg import pinv,norm
from snowdrop.src.utils.prettyTable import PrettyTable
NITERATIONS = 100 # Maximum nymber of iterations
EPSILON = 1.e-5 # Convergence tolerance
Nt = 10
Ngrid = 100
ngrid = 10
# Parameters
period_duration = 60/2 # We are assuming that average life expectancy is 60 years
alpha = 0.5 # Cobb-Douglas productivity function parameter, exponent
beta = 0.5 # Impatience parameter
pi = 10 # Cobb-Douglas productivity function parameter, normalization
phi = 0.9 # Selfishness parameter
labor = 1 # Unskilled labor multiplier
U = 1.03**period_duration # Final machine productivity
u = np.ones(Nt); u[3:] = U # Machine productivity
R = alpha*pi*((1-alpha)/alpha/U)**(alpha-1)
# Initialize variables
L = 1
N=np.zeros(Nt); Q=np.ones(Nt); QL=np.zeros(Nt); QM=np.zeros(Nt); QS=np.zeros(Nt); WL=np.zeros(Nt)
I1=np.zeros(Nt); I2=np.zeros(Nt); C1=np.zeros(Nt); C2=np.zeros(Nt); Sv=np.zeros(Nt)
B=np.zeros(Nt); U1=np.zeros(Nt); U2=np.zeros(Nt); M=np.ones(Nt); S=np.ones(Nt)
[docs]
def value_iteration(Output=False):
"""
Solve Overlapping Generations Model by value function iteration.
Parameters:
Output: If True save results in excel file.
Returns:
iterations : int
Number of iterations.
crit : float
Convergence value.
y : numpy array
Solution of Bellman equation
"""
global u, M, S
crit = 1.e6 # convergence value
iterations = 0 # number of iterations
dr = np.zeros((Ngrid,Ngrid),dtype=int) # decision rule
# Initial guess for value function
Tv = np.zeros((Ngrid,Ngrid))
utility = np.zeros((Ngrid,Ngrid))
x = np.zeros((Ngrid,Ngrid))
v = 0
sv = 0
# Need to seed M, typically with 'equilibrium' values.
Mt,St = MSequilibrium(L,u[0]);
M *= Mt; S *= St
# Define grid
cy = np.linspace(1.e-6,10,Ngrid)
co = np.linspace(1.e-6,100,Ngrid)
for t in range(1,Nt):
# The marginal products of current M,S determine the income/consumption
# of the current old generation.
n = L + u[t]*Mt
s = (1-alpha)*(sv+L/u[t])
QM = dQdM(n,s,u[t])
QS = dQdS(n,s)
### Solve Bellamn equation
while crit > EPSILON:
m = Mt; s = St; sv = 0
# Utility
for i in range(Ngrid):
cy_g = cy[i]
for j in range(Ngrid):
co_g = co[j]
utility = phi*beta*np.log(cy_g) + phi*(1-beta)*np.log(co_g)
x[i,j] = utility + (1-phi)*v
# find new policy rule
Tv[i,j] = np.max(x)
dr[i,j] = np.argmax(x)
# Decision rules
for i in range(Ngrid):
for j in range(Ngrid):
cy = cy[dr[i,j]]
co = co[dr[i,j]]
# Update variables
n = L + u*m
# Skills
s = (1-alpha)*(sv+L/u[t])
# Wage rate
w = alpha * pi * (n/s)**(alpha-1)
# Income of old
io = QM*(m+s)
# Bequest from old to young
b = np.max(0,io - co)
# Income of young
iy = w*L + b
# Savings of youn
sv = iy - cy
# sv = m + s
# Perfect foresigh assumption
m = alpha*sv + (1-alpha)*L/u
s = (1-alpha)/alpha*n
crit= np.max(abs(Tv-v))/norm(Tv)
v = Tv.copy()
iterations += 1
print('Iteration # {} \tCriterion: {:.2e}'.format(iterations,crit))
# Compute intermediate product N:
N[t] = L + u[t]*M[t]
# Compute productivity and marginal products:
Q[t] = Qtot(N[t],S[t])
QL[t] = dQdL(N[t],S[t])
WL[t] = QL[t]*L
# Compute interest factor. Use market clearing condition:
R = QM[t]
# Compute constant m:
m = WL[t]*R/(R-1)
# Old generation income is a return on machines and skilled labor
I2[t] = QM[t]*M[t] + QS[t]*S[t]
# Compute bequest of old generation.
bequest = (1-phi)/(1-phi*beta)*I2[t] - m*phi*(1-beta)/(1-phi*beta)
# It can not become negative.
b[t] = min(I2[t],max(0,bequest))
# Compute consumption of old generation
C2[t] = I2[t] - b[t]
# Consumption of young generation:
# Income of young generation is labor plus bequest
I1[t] = WL[t] + b[t]
C1[t] = beta*phi*(I1[t] + m/R)
# Consumption can not be greater than income
C1[t] = min(I1[t],C1[t])
# Savings of young generation:
Sv[t] = I1[t] - C1[t]
# Perfect foresight assumption means the Young now guess the next M and
# S values, using their unanticipated guess for u. These become the
# values of M,S for next period.
Mt,St = MSnext(L[t],Sv[t],u[t])
M[t+1] = max(0,Mt)
S[t+1] = max(0,St)
# Utility definition
U1[t] = beta*phi*np.log(C1[t]) # This is the current young utility
U2[t] = (1-beta)*phi*np.log(C2[t]) # This is the current old utility
# Compute Lifetime Utility as a sum of utilities for each generation:
LU = list(U1[:-1] + U2[1:])
LU.append(0)
U = np.exp(LU)
share_of_labor = (QL*L+QS*S)/Q
# Display results:
table = [m, S, QL, QS, C1, C2, Sv, Q, I1, I2, share_of_labor, LU, U, b]
table = np.around(table, decimals=2)
table = np.row_stack((u, table))
table = np.row_stack((np.arange(Nt), table))
header = ['T','u','m','Mg','Mp','TP','S','Wl','Ws','C','D','Sv','Q','Iy','Io','Share of Labor','LU','U','Bequest']
pTable = PrettyTable(header)
pTable.float_format = '6.2'
for i in range(1,Nt-2):
pTable.add_row(table[:,i])
if Output:
print()
print('Bequest model:')
print(pTable)
print()
table = table.T
path = os.path.dirname(os.path.abspath(__file__))
fname = path + '\\..\\data\\bequest.csv'
rows,columns = table.shape
with open(fname, 'w') as f:
f.writelines(",".join(header) + "\n")
for r in range(1,rows-2):
for c in range(columns):
if c == columns:
f.writelines(table[r,c])
else:
f.write(str(table[r,c]) + ",")
f.write("\n")
var_names
var_labels
return iterations,crit,y,var_names,var_labels
[docs]
def MSequilibrium(Labor,u):
# Find equilibrium M,S values.
A = np.array([[1, 1], [-1, alpha/(1-alpha)]])
WL = pi*alpha*(alpha*u/(1-alpha))**(alpha-1)
B = np.array([(1-beta)*WL, Labor/u])
M,S = np.linalg.solve(A,B)
return M,S
[docs]
def MSnext(Labor,Saving,u):
# Find non-equilibrium M,S values.
A = np.array([[1, 1], [-1, alpha/(1-alpha)]])
B = np.array([ Saving, Labor/u ])
M,S = np.linalg.solve(A,B)
return M,S
[docs]
def Qtot(N,S):
# Compute total productivity Q.
Q = pi * N**alpha * S**(1-alpha)
return Q
[docs]
def dQdL(N,S):
# Compute marginal product w.r.t. L.
QL = alpha * pi * (N/S)**(alpha-1)
return QL
[docs]
def dQdM(N,S,u):
# Compute marginal product w.r.t. M.
QM = alpha * pi * u * (N/S)**(alpha -1)
return QM
[docs]
def dQdS(N,S):
# Compute marginal product w.r.t. S.
QS = (1-alpha) * pi * (N/S)**(alpha )
return QS
[docs]
def Plot(path_to_dir,data,variable_names,var_labels={}):
from snowdrop.src.graphs.util import plotTimeSeries
header = 'Solution of Bellman Equation'
series = []; labels = []; titles = []
if np.ndim(data) == 3:
ns,n,nvar = data.shape
for i in range(nvar):
ser = []; lbls = []
for j in range(ns):
ser.append(pd.Series(data[j,:,i]))
lbls.append("Markov State " + str(1+j))
series.append(ser)
labels.append(lbls)
var = variable_names[i]
name = var_labels[var] if var in var_labels else var
titles.append(name)
elif np.ndim(data) == 2:
n,nvar = data.shape
for i in range(nvar):
ser = pd.Series(data[:,i])
series.append([ser])
labels.append([])
var = variable_names[i]
name = var_labels[var] if var in var_labels else var
titles.append(name)
else:
nvar = len(y)
for i in range(nvar):
ser = []; lbls = []
data = y[i]
ns = len(data)
for j in range(ns):
ser.append(pd.Series(data[j]))
if ns>1:
lbls.append("Markov State " + str(1+j))
else:
lbls.append(" ")
series.append(ser)
labels.append(lbls)
var = variable_names[i]
name = var_labels[var] if var in var_labels else var
titles.append(name)
if nvar <= 4:
sizes = [2,2]
else:
sizes = [int(np.ceil(nvar/3)),3]
plotTimeSeries(path_to_dir=path_to_dir,header=header,titles=titles,labels=labels,series=series,sizes=sizes)
if __name__ == '__main__':
"""The main entry point."""
# Solve Bellman equation
t0 = time()
iterations,crit,y,var_names,var_labels = value_iteration()
elapsed = time() - t0
print("\nElapsed time: %.2f (seconds)" % elapsed)
Plot(path_to_dir=path_to_dir,data=y,variable_names=variable_names,var_labels=var_labels)