Problem Set 3: Problem 3 Coding Solutions

using LinearAlgebra
using Plots
using Pkg
#Pkg.add("Distributions")
using Distributions
using Random
using Statistics

Part 1

Part (a) : policy function for an infinitely lived agent

Using solutions from problem set 2 problem 4

## Specify parameters values ##
δ = 0.9

# Make vector with \tilde{y} values and probabilities (Py)
Y = [y for y in 0.1:0.1:0.9]
Py = ones(9)./9 # uniform distribution probabilities

# Wealth level vector #
X = [x for x in 0.01:0.01:10];

## Specify initial guess ##
V0 = X;

#### Create two matricies for each set of pairs (s, y) \in (X x Y) ####
# Note that this gives us two 1000x9 matrities #
s_vals = repeat(X,1,9)
y_vals = repeat(Y',1000,1);

# Calculate x_{+1} = s + \tilde{y}. Note, however, that we don't allow for x_{+1} > 10 because that is outside our pre-defined grid #
new_vals = min.(s_vals + y_vals,10);

## For a given value function, calculate VC ##
function VC_function(V0)
    V_vals = V0[Int.(round.(new_vals*100))] 
    VC = V_vals*Py 
    return VC
end

X_vals_2 = repeat(X,1,1000)
s_vals = X_vals_2';

c_bar = 1e-10 # set c_bar given in problem

# Gives us a 1000 x 1000 matrix for a given par (x,s) #
U_vals = log.(max.(X_vals_2-s_vals,c_bar));

function iterate_V(V;policy = false)
    # For a given initial value function V, calculate the continuation value using VC_function defined above #
    # Input: Value function V over every possible value of s i.e. it is a 1000x1 vector #
    # Output w/ repeat function: 1000 x 1000 matrix with each column corrresponding to a value of (s) 
    VC_vals = repeat(VC_function(V)',1000,1) 

    # Add flow utility + continuation value #
    V1_vals = U_vals + δ*VC_vals

    # Over every wealth of x (i.e. across every row since each row corresponds with a value of x), find the max value function  
    # value V1 AND also find the index location of where the max value occurs (this is important for using being able to find the amount 
    # of savings s that maximizes the value function)
    V1, temp = findmax(V1_vals,dims=2)
    if policy == true
        # Remember the relationship between the column index and the value of saving $s$. Here, we extract the column index for each row where 
        # the max value function occurs. And then, we divide by 100 which gives us the value of saving $s$ that maximizes the value function.
        return [s[2] for s in temp]./100 
    end
    return V1
end

diff = 1000 # intialize difference variable
tol = 1e-8 # set tolerance as defined in problem

## Keep iterating until the difference between value function iterations gets very small ##
while diff > tol
    V1 = iterate_V(V0)
    diff = norm(V1-V0)
    V0=V1
end

V_inf = V0;
s_inf = iterate_V(V_inf, policy = true);
# Manually set saving when x = 0.01 to 0 - otherwise negative consumption #
s_inf[1] = 0;

Part (b) : policy functions of an agent who faces the same problem but lives only for three periods

There are three possible ages, so there are three value functions. This also means that we will have three policy functions.

Part (c): using policy functions that we derived in problem 1 of this problem set

In this problem, $T = 3$. Note that we derived that $$ c_{T-1} = \frac{1}{1+\delta} x $$

$$ c_{T-2} = \frac{1}{1+\delta + \delta^2} x $$

We know that in the terminal period in this model, the agent will consume everything so $c_T = x$.

Also note that the policy function we are using is saving so $s = x - c$

#### Part (c) : define X1, X2, and X3 ####
# These give us the saving policy rule for each age = 1, = 2, = 3 #

X1 = X - (1/(1+δ+δ^2))*X; # age = 1 (or t = T -2)
X2 = X - (1/(1+δ))*X; # age = 2 (or t = T -1)
X3 = X - X; # age = 3 (or t = T)

s_finite = [X1 X2 X3];

Part 2: Implications for consumption following unexpected gov transfer

Think about MPCs of these policies.

plot(X, X1, label="X1 (age=1)")
plot!(X, X2, label="X2 (age=2)")
plot!(X, X3, label="X3 (age=3)")

Part 3: Equilibrium distribution of wealth

For this question, we are asked to find the equilibrium distirubtion of wealth for the infinite-lived agent model and for the finite-lived agent model.

Infinetly-Lived Agent

#### Subpart (a) : simulate 12,000 agents from uniform distribution of wealth ####

N = 12000 # Number of agents
W0 = rand(X,N) # Initial distribution of wealth, we assume it's uniform over X
W0_initial = W0 # saving for plot

#### Subpart (b) : Draw income shock from Y ####

Y = [y for y in 0.1:0.1:0.9] 
Y_draws = rand(Y,N)

#### Subpart (c) & (d) : Function that takes wealth and returns optimal savings S ####

function s(w)
    w_index = Int(round(w*100)) # For a given amount of wealth w
    # Find the optimal saving from the policy function (remember the trick about relationship between index and X value)
    savings = s_inf[min(w_index,1000)] 
    return savings
end

#### Subpart (e) & (f) : Repeat process for savings over time for T = 10000  ####
# For this problem, it is as if we have 10,000 time periods where in each period, the agent 
# gets some draw of y and savings a certain amount. 

T = 10000 # number for periods
for t in 1:T # for each period
    yt = rand(Y,N) # recieve a draw of y
    # the wealth in the next period for that agent is the amount they save s(W0[i]) + their draw of income
    W1 = [s(W0[i]) + yt[i] for i in 1:N] 
    W0 = W1 # W1 is the new distribution of wealth
end

W_inf = W0 ;

#### Subpart (g) : Plot distributions ####

histogram(W_inf, label="W_inf", alpha=0.5, title="Wealth Distribution", xlabel="Wealth", ylabel="Frequency")
#histogram!(W0_initial, label="W0_initial", alpha=0.5, bins=30)

Finite-Lived Agent

Very similar code as above but now have savings indexed by age

#### Subpart (a) : simulate 12,000 agents from uniform distribution of wealth ####
N = 12000 #Number of agents

#Age vector
age = ones(N) 
age[4001:8000] .= 2 
age[8001:12000] .= 3 
age = Int.(age)

W0 = rand(X,N) #randomly drawing wealth from X

#### Subpart (a) : simulate 12,000 agents from uniform distribution of wealth ####
W0_fin = [W0 age]

#### Subpart (c) & (d) : Function that takes wealth and returns optimal savings S ####

function s_fin(w,a)
    w_index = Int(round(w*100))
    # Find the optimal saving from the policy function (remember the trick about relationship between index and X value)
    savings = s_finite[min(w_index,1000),Int(a)]
    return [savings a] 
end 

#### Subpart (e) & (f) : Repeat process for savings over time for T = 10000  ####
# For this problem, it is as if we have 10,000 time periods where in each period, the agent 
# gets some draw of y and savings a certain amount. 

T = 10000 # number for periods

for t in 1:T # for each period
    yt = rand(Y,N) # recieve a draw of y
    W1 = zeros(N, 2)  # Initialize W1 as a matrix to store wealth and age
    
    for i in 1:N 
        # the wealth in the next period for that agent - note first argument of s_fin is wealth amount and second is age #
        result = s_fin(W0_fin[i, 1], W0_fin[i, 2])  # This returns a 1x2 matrix [savings, age]
        W1[i, 1] = result[1] + yt[i]  # New wealth for the next period
        W1[i, 2] = result[2] # age in period t

        # Increase age by one. If age = 3, then replace with youngest agent #
        if  W1[i, 2] <= 2
                W1[i, 2] = W1[i, 2] + 1 
        else # replace with youngest agent and draw initial wealth randomly 
            W1[i, 2] = 1
            W1[i, 1] = rand(W_inf)
        end
    end 
    
    W0_fin = W1 # W1 is the new distribution of wealth
end

W_fin = W0_fin ;

#### Subpart (g) : Plot distributions ####

histogram(W_fin[:,1], label="W_finite", alpha=0.5, title="Wealth Distribution", xlabel="Wealth", ylabel="Frequency")
#histogram!(W0_initial, label="W0_initial", alpha=0.5)

### Comparing distribution of finite and infinte agents ####

histogram(W_fin[:,1], label="W_fin", alpha=0.5, title="Wealth Distribution", xlabel="Wealth", ylabel="Frequency")
histogram!(W_inf, label="W_inf", alpha=0.5)

Part 4: Compute the change in consumption

The change in consumption is just $\Delta C_i = C_i (W + 1) - C_i (W)$.

Note that we can calculate the aggregate MPC as

$$ MPC = \frac{\sum \Delta C_i}{1 \cdot R} $$

where $R$ is the number of recipents of the stimulus.

### Function to calculate the change in consumption ###

function c_change(w;inf = true,age = 0)
    if inf == true
        original_c = w - s(w)
        new_c = w + 1 - s(w+1)
    elseif inf == false
        original_c = w - s_fin(w,age)[1]
        new_c = w + 1 - s_fin(w+1,age)[1]
    end
    return new_c - original_c
end

c_change (generic function with 1 method)

ΔC_inf = zeros(N)
ΔC_fin = zeros(N)

## Change in consumption for the infinite-lived model ##
for i in 1:N
    if W_inf[i] < median(W_inf) 
        ΔC_inf[i] = c_change(W_inf[i]) # every agent below the median gets a stimulus 
    end
end

## Change in consumption for the finite-lived model ##
for i in 1:N
    if W_fin[i, 1] < median(W_fin[:, 1]) 
        ΔC_fin[i] = c_change(W_fin[i, 1], inf = false, age = W_fin[i, 2]) # every agent below the median gets a stimulus 
    end
end

## Plot distribution of ΔC ##
histogram(ΔC_inf,alpha=0.5, xlabel="ΔC", ylabel="Frequency")
histogram!(ΔC_fin,alpha=0.5)

## Calculate the aggregate MPC ##

# Infinite-lived agents #
agg_spend = sum(W_inf[i] < median(W_inf)  for i in 1:N)
mpc_inf = sum(ΔC_inf)/agg_spend

# Finite-lived agents #
agg_spend = sum(W_fin[i, 1] < median(W_fin[:, 1])  for i in 1:N)
mpc_fin = sum(ΔC_fin)/agg_spend

println("mpc_inf: ", mpc_inf)
println("mpc_fin: ", mpc_fin)