Guvenen (2009): Asset Pricing with Heterogeneous IES¶
Guvenen (2009) constructs a two-agent model to account for several salient features of asset pricing moments, such as high risk premium, low and relatively smooth interest rate, and countercyclical movements in risk premium and Sharpe ratio. The two key ingredients of his model are limited stock market participation and heterogeneity in the elasticity of inter-temporal substitution in consumption (EIS).
When all endogenous states are pre-determined, such as capital and bond share in the current example, the expectation of future variables can be parameterized and integrated before solving the system of equilibrium equations, which greatly speeds up the computation. We illustrate how to use parameterized expectation with the toolbox in this example.
The Model¶
There are two types of infinitely-lived agents: stockholders (\(h\)) with measure \(\mu\), and non-stockholders (\(n\)) with measure \(1-\mu\). Agents have Epstein-Zin utility functions
for \(i=h,n\). Most importantly, \(\rho^{h}<\rho^{n}\), i.e., the non-stockholders have lower EIS which is inversely proportional to \(\rho^{i}\), and thus they have higher desire for consumption smoothness over time. Each agent has one unit of labor endowment.
Stockholders can trade stock \(s_{t}\) and bond \(b_{h,t}\) at prices \(P_{t}^{s}\) and \(P_{t}^{f}\) respectively. Their budget constraint is
where \(W_t\) is the labor income and borrowing constraint is
and in calibration \(\underline{B}\) is set to six times of the average monthly wage rate. The non-stockholders have the same constraints. In addition, they are restricted from trading stocks.
A representative firm produces the consumption good using capital \(K_{t}\) and labor \(L_{t}\) based on a Cobb-Douglas production function:
and the technology evolves according to an AR(1) process:
The firm maximizes its value \(P_{t}^{s}\) expressed as the sum of its future dividends \(\left\{ D_{t+j}\right\} _{j=1}^{\infty}\) discounted by the shareholders’ marginal rate of substitution process:
The firm accumulates capital subject to a concave adjustment cost function in investment:
Each period, the firm sells one-period bonds at price \(P_{t}^{f}\). The bond supply is constant and equals to \(\chi\) fraction of its average capital stock \(\bar{K}\). Thus dividend \(D_{t}\) can be written as
A sequential competitive equilibrium is given by sequences of allocations \(\{c_{i,t},b_{i,t+1},s_{t+1},I_{t},K_{t+1},L_{t}\}\) \(i=h,n\) and prices \(\left\{ P_{t}^{s},P_{t}^{f},W_{t}\right\}\) such that (i) given the price sequences, \(\left\{ c_{i,t},b_{i,t+1},s_{t+1}\right\}\) \(i=h,n\) solve the stockholders’ and non-stockholders’ optimization problems; (ii) Given the wage sequence \(\left\{ W_{t}\right\}\) and the law of motion for capital, \(\left\{L_{t},I_{t}\right\}\) are optimal for the representative firm; (iii) all markets clear:
We compute a recursive equilibrium using \(\left\{ K_{t},B_{t}^{n},Z_{t}\right\}\) as the aggregate state variables, where \(B_{t}^{n}=\left(1-\mu\right)b_{n,t}\) is total bond holding by the non-stockholders. We have \(8\) unknowns: {\(c_{h,t}\), \(c_{n,t}\), \(I_{t}\), \(B_{t+1}^{n}\), \(\lambda_{h,t}\), \(\lambda_{n,t}\), \(P_{t}^{s}\), \(P_{t}^{f}\)}, and the \(8\) equations used to solve them are:
Euler equations for bond holding:
Euler equations for the stockholders’ demand of equity:
Slackness condition of borrowing limit:
The budget constraints (imposing \(s_{t+1}=1/\mu\)):
Firm’s optimal capital accumulation \(K_{t+1}\):
in which capital price \(q_{t}\) is the Lagrangian multiplier on the capital formation equation and satisfies
The auxiliary variables can be determined by the utility function, market clearing conditions, and the following two equations:
In period \(t\), the 6 future variables in use: \(c_{h,t+1}\), \(c_{n,t+1}\), \(P_{t+1}^{s}+D_{t+1}\), \(I_{t+1}/K_{t+1}\), \(U_{h,t+1}\) and \(U_{n,t+1}\) are functions of \(\left\{ K_{t+1},B_{t+1}^{n},Z_{t+1}\right\}\) and are solved from the previous iteration. Similar to Guvenen (2009), the initial guesses for these functions are obtained by solving a version of the model with no leverage.
The gmod File¶
The model can be solved with the following guvenen2009.gmod
file
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 | % Parameters
parameters beta alpha rhoh rhon theta delta mu xsi chi a1 a2 Kss Bbar bn_shr_lb bn_shr_ub varianceScale;
beta = 0.9966; % discount factor
alpha = 6; % risk aversion
rhoh = 1/.3; % inv IES for stockholders
rhon = 1/.1; % inv IES for non-stockholders
theta = .3; % capital share
delta = .0066; % depreciation rate
mu = .2; % participation rate
xsi = .4; % adjustment cost coefficient
chi = .005; % leverage ratio
a1 = (((delta^(1/xsi))*xsi)/(xsi-1));
a2 = (delta/(1-xsi));
Kss = ((1/beta-1+delta)/theta)^(1/(theta-1));
Bbar = -0.1*(1-theta)*Kss^theta; %borrowing constraint
varianceScale = 1e4;
TolEq = 1e-4;
INTERP_ORDER = 4; EXTRAP_ORDER = 4;
PrintFreq = 100;
SaveFreq = inf;
% Shocks
var_shock Z;
shock_num = 15;
phi_z = 0.984; % productivity AR(1)
mu_z = 0;
sigma_e = 0.015/(1+phi_z^2+phi_z^4).^0.5;
[z,shock_trans,~]=tauchen(shock_num,mu_z,phi_z,sigma_e,2);
Z = exp(z);
% States
var_state K bn_shr;
K_pts = 10;
K = exp(linspace(log(.84*Kss),log(1.2*Kss),K_pts));
bn_shr_lb = (1-mu)*Bbar/(chi*Kss);
bn_shr_ub = (chi*Kss - mu*Bbar)/(chi*Kss);
b_pts = 30;
bn_shr = linspace(bn_shr_lb,bn_shr_ub,b_pts);
% Last period
var_policy_init c_h c_n;
inbound_init c_h 1e-6 100;
inbound_init c_n 1e-6 100;
var_aux_init Y W vh vn vhpow vnpow Ps Pf Div Eulerstock Eulerbondh Eulerbondn Inv dIdK Eulerf;
model_init;
Y = Z*(K^theta);
W = (1-theta)*Z*(K^theta);
resid1 = 1 - (W + (bn_shr*chi*Kss/(1-mu)))/c_n; % c_n: individual consumption
resid2 = 1 - (W + (Div/mu) + ((1-bn_shr)*chi*Kss/mu))/c_h; % c_h: individual consumption
vh = ((1-beta)*(c_h^(1-rhoh)))^(1/(1-rhoh));
vn = ((1-beta)*(c_n^(1-rhon)))^(1/(1-rhon));
vhpow = vh^(1-alpha);
vnpow = vn^(1-alpha);
Pf = 0;
Ps = 0;
Div = Y - W - (1-Pf)*chi*Kss; % investment is zero
Eulerstock = (vh^(rhoh-alpha))*(c_h^-rhoh)*(Ps + Div);
Eulerbondh = (vh^(rhoh-alpha))*(c_h^-rhoh);
Eulerbondn = (vn^(rhon-alpha))*(c_n^-rhon);
Inv = 0;
Knext = 0;
dIdK = (Inv/K) - (1/a1)*(xsi/(xsi-1))*(Inv/(K*((1/a1)*((Knext/K)-(1-delta)-a2))))*(Knext/K);
Eulerf = (vh^(rhoh-alpha))*(c_h^-rhoh)*(theta*Z*(K^(theta-1)) - dIdK);
equations;
resid1;
resid2;
end;
end;
var_interp EEulerstock_interp EEulerbondh_interp EEulerbondn_interp EEulerf_interp Evh_interp Evn_interp EPD_interp EPD_square_interp;
initial EEulerstock_interp shock_trans*reshape(Eulerstock,shock_num,[]);
initial EEulerbondh_interp shock_trans*reshape(Eulerbondh,shock_num,[]);
initial EEulerbondn_interp shock_trans*reshape(Eulerbondn,shock_num,[]);
initial EEulerf_interp shock_trans*reshape(Eulerf,shock_num,[]);
initial Evh_interp shock_trans*reshape(vhpow,shock_num,[]);
initial Evn_interp shock_trans*reshape(vnpow,shock_num,[]);
initial EPD_interp shock_trans*reshape(Div,shock_num,[]);
initial EPD_square_interp shock_trans*reshape(Div.^2,shock_num,[]) / varianceScale;
EEulerstock_interp = shock_trans*Eulerstock;
EEulerbondh_interp = shock_trans*Eulerbondh;
EEulerbondn_interp = shock_trans*Eulerbondn;
EEulerf_interp = shock_trans*Eulerf;
Evh_interp = shock_trans*vhpow;
Evn_interp = shock_trans*vnpow;
EPD_interp = shock_trans*(Ps+Div);
EPD_square_interp = shock_trans*(Ps+Div).^2 / varianceScale;
% Endogenous variables, bounds, and initial values
var_policy c_h c_n Ps Pf Inv bn_shr_next lambdah lambdan;
inbound c_h 1e-3 10;
inbound c_n 1e-3 10;
inbound Ps 1e-3 500;
inbound Pf 1e-3 10;
inbound Inv 1e-9 10;
inbound bn_shr_next bn_shr_lb bn_shr_ub;
inbound lambdah 0 2;
inbound lambdan 0 2;
% Other equilibrium variables
var_aux Y W b_h b_n Div dIdKp Eulerstock Eulerbondh Eulerbondn dIdK Eulerf vhpow vnpow omega PDratio Rs R_ep vh vn Knext std_ExcessR SharpeRatio;
model;
Y = Z*(K^theta); % output
W = (1-theta)*Z*(K^theta); % Wage = F_l
Div = Y - W - Inv - (1-Pf)*chi*Kss; % dividends
Knext = (1-delta)*K + (a1*((Inv/K)^((xsi-1)/xsi))+a2)*K;
dIdKp = (1/a1)*(xsi/(xsi-1))*(Inv/(K*((1/a1)*((Knext/K)-(1-delta)-a2))));
b_h = (1-bn_shr)*chi*Kss/mu;
b_n = bn_shr*chi*Kss/(1-mu);
[EEulerstock_future,EEulerbondh_future,EEulerbondn_future,EEulerf_future,Evh_future,Evn_future,EPD_future,EPD_square_future] = GDSGE_INTERP_VEC(shock,Knext,bn_shr_next);
EPD_square_future = EPD_square_future*varianceScale;
vh = ((1-beta)*(c_h^(1-rhoh)) + beta*(Evh_future^((1-rhoh)/(1-alpha))))^(1/(1-rhoh));
vn = ((1-beta)*(c_n^(1-rhon)) + beta*(Evn_future^((1-rhon)/(1-alpha))))^(1/(1-rhon));
Eulerstock = (vh^(rhoh-alpha))*(c_h^-rhoh)*(Ps + Div);
Eulerbondh = (vh^(rhoh-alpha))*(c_h^-rhoh);
Eulerbondn = (vn^(rhon-alpha))*(c_n^-rhon);
dIdK = (Inv/K) - (1/a1)*(xsi/(xsi-1))*(Inv/(K*((1/a1)*((Knext/K)-(1-delta)-a2))))*(Knext/K);
Eulerf = (vh^(rhoh-alpha))*(c_h^-rhoh)*(theta*Z*(K^(theta-1)) - dIdK);
vhpow = vh^(1-alpha);
vnpow = vn^(1-alpha);
omega = (Ps+Div+ mu*b_h)/(Ps+Div+chi*Kss);
PDratio = Ps/Div;
Rs = EPD_future/Ps;
R_ep = Rs - 1/Pf;
% The following inline implements
% std_ExcessR = (GDSGE_EXPECT{(PD_future'/Ps - Rs)^2})^0.5;
std_ExcessR = (EPD_square_future/(Ps^2) + Rs^2 - 2*EPD_future*Rs/Ps)^0.5;
SharpeRatio = R_ep/std_ExcessR;
% Equations:
err_bdgt_h = 1 - (W + (Div/mu) + b_h - Pf*(chi*Kss*(1-bn_shr_next)/mu))/c_h; % these are individual consumptions
err_bdgt_n = 1 - (W + b_n - Pf*(bn_shr_next*chi*Kss/(1-mu)))/c_n;
foc_stock = 1 - (beta*EEulerstock_future*(Evh_future^((alpha-rhoh)/(1-alpha))))/((c_h^(-rhoh))*Ps);
foc_bondh = 1 - (beta*EEulerbondh_future*(Evh_future^((alpha-rhoh)/(1-alpha))) + lambdah)/((c_h^(-rhoh))*Pf);
foc_bondn = 1 - (beta*EEulerbondn_future*(Evn_future^((alpha-rhon)/(1-alpha))) + lambdan)/((c_n^-rhon)*Pf);
foc_f = 1 - (beta*EEulerf_future*(Evh_future^((alpha-rhoh)/(1-alpha))))/((c_h^(-rhoh))*dIdKp);
slack_bn = lambdan*(bn_shr_next - bn_shr_lb); %mun_lw*bn_shr_next;
slack_bh = lambdah*(bn_shr_ub - bn_shr_next); %mun_up*(1-bn_shr_next);
equations;
err_bdgt_h;
err_bdgt_n;
foc_stock;
foc_bondh;
foc_bondn;
foc_f;
slack_bn;
slack_bh;
end;
end;
simulate;
num_periods = 10000;
num_samples = 100;
initial K Kss;
initial bn_shr 0.5;
initial shock 2;
var_simu Y c_h c_n Inv Ps Div Pf bn_shr_next Knext omega PDratio Rs R_ep SharpeRatio std_ExcessR;
K' = Knext;
bn_shr' = bn_shr_next;
end;
|
The implementation uses the exact parameterization and specification (number of grid points, domains of endogenous state variables, number of discretized exogenous shocks, the cubic spline approximation, etc.) in Guvenen (2009). However, the algorithm differs from his. His is based on the algorithm in Krusell and Smith (1998): starting from a conjectured law of motion for state-variables and pricing functions, he solves the agents’ Bellman equation and the agents’ policy functions using standard value function iterations. Then he uses these policy functions and temporary market clearing conditions to obtain a new law of motions and new pricing functions. The current algorithm, like other examples, utilizes the fact that agents’ optimization is concave, replaces the agent’s optimization with first order and complementarity conditions, and solves the agents’ optimization problem together with the market clearing conditions.
Since both state variables \(K_t\) and \(B_t^n\) are predetermined, the expectation of future variables (e.g., the marginal utility entering the bond’s Euler equation, or the price and dividend entering the equity’s Euler equation) can be calculated before-hand conditional on current exogenous states. For example, for the bond’s Euler equation
we can define function
and approximate \(\tilde{\Lambda}(z,K',B')\) instead of \(\Lambda(z',K',B')\). This is made possible by that \((K,B)\) are predetermined in that the transition of the two states do not depend on \(z'\) (this is not true when the endogenous state follows an implicit state transition function that depends on the realization of \(z'\), such as wealth share and consumption share). When a projection method or spline is used to approximate \(\tilde{\Lambda}(z,K,B)\), the integration becomes quite simple, since projection is a linear operator and therefore the projection of an expectation equals to the expectation of the projection.
Once \(\tilde{\Lambda}(z,K',B')\) is approximated beforehand, then one only needs to evaluate \(\tilde{\Lambda}\) once for the current exogenous shock (\(z\)) in each evaluation of the system of equations. If instead one is to work with \(\Lambda(z',K',B')\), \(\Lambda\) needs to be evaluated for each \(z'\) and integrated using the transition matrix. This is why integrating \(\Lambda(z',K',B')\) beforehand is more computation-efficient, especially when the size of exogenous shocks is large.
In particular, the parameterized expectation is implemented by declaring the expectation of relevant variables in var_interp in
79 | var_interp EEulerstock_interp EEulerbondh_interp EEulerbondn_interp EEulerf_interp Evh_interp Evn_interp EPD_interp EPD_square_interp;
|
specifying the updating of var_interp by integrating the corresponding variables using the transition matrix
89 90 91 92 93 94 95 96 | EEulerstock_interp = shock_trans*Eulerstock;
EEulerbondh_interp = shock_trans*Eulerbondh;
EEulerbondn_interp = shock_trans*Eulerbondn;
EEulerf_interp = shock_trans*Eulerf;
Evh_interp = shock_trans*vhpow;
Evn_interp = shock_trans*vnpow;
EPD_interp = shock_trans*(Ps+Div);
EPD_square_interp = shock_trans*(Ps+Div).^2 / varianceScale;
|
and evaluating the approximation at the current shock only when evaluating the system of equations
124 | [EEulerstock_future,EEulerbondh_future,EEulerbondn_future,EEulerf_future,Evh_future,Evn_future,EPD_future,EPD_square_future] = GDSGE_INTERP_VEC(shock,Knext,bn_shr_next);
|
where the GDSGE_INTERP_VEC takes the first argument shock, a reserved keyword referring to the index of the current exogenous shock.
After processed by the online compiler, we can run the returned file to run policy iterations:
>> options = struct;
options.MaxIter = 500;
IterRslt = iter_guvenen2009(options);
Iter:100, Metric:0.512539, maxF:8.65033e-09
Elapsed time is 8.272149 seconds.
Iter:200, Metric:0.285722, maxF:9.39891e-09
Elapsed time is 7.549686 seconds.
Iter:300, Metric:0.174423, maxF:9.01752e-09
Elapsed time is 7.250628 seconds.
Iter:400, Metric:0.113139, maxF:8.89712e-09
Elapsed time is 6.938417 seconds.
Iter:500, Metric:0.0763611, maxF:8.61954e-09
Elapsed time is 6.762700 seconds.
We first solve a simple version of the problem as a warm-up. This simple problem has a lower aggregate net supply of bond (\(\chi=0.005\)) and a tighter borrowing constraint (\(\bar{B}\)). This ensures that starting from a last period solution the equilibrium exists, which is also noted in Guvenen (2009).
Starting with the warm-up solution, we now change the parameters, increase \(\chi\) to \(0.15\), its value in the paper, and specify a much more relaxed borrowing constraint (\(\bar{B}\)). This can be done very conveniently with the toolbox by specifying the new parameters into a structure, and the warm-up solution as a field named WarmUp into the structure.
>> beta = 0.9967; % discount factor
theta = .3; % capital share
delta = .0066; % depreciation rate
mu = .2; % participation rate
Kss = ((1/beta-1+delta)/theta)^(1/(theta-1));
chi = 0.15;
Bbar = -6*(1-theta)*Kss^theta;
bn_shr_lb = (1-mu)*Bbar/(chi*Kss);
bn_shr_ub = (chi*Kss - mu*Bbar)/(chi*Kss);
b_pts = 30;
options = struct;
options.chi = chi;
options.Bbar = Bbar; %borrowing constraint
options.bn_shr_lb = bn_shr_lb;
options.bn_shr_ub = bn_shr_ub;
options.bn_shr = linspace(bn_shr_lb,bn_shr_ub,b_pts);
options.WarmUp = IterRslt;
options.SkipModelInit=1;
We can now call the routine with the option structure as an argument
>> IterRslt = iter_guvenen2009(options);
Iter:600, Metric:0.0605821, maxF:9.85466e-09
Elapsed time is 6.361178 seconds.
...
Iter:2490, Metric:9.99635e-05, maxF:9.29708e-09
Elapsed time is 5.205234 seconds.
The policy iteration continues from where the warm-up solution ends, toward convergence.
We can inspect the policy functions at convergence by calling in MATLAB:
figure; clf;
subplot(1,2,1);
surf(KMesh,BMesh,squeeze(var_aux.R_ep(z1,:,:))*12);
title('Equity Premium');
xlabel('$K$','Interpreter','Latex');
ylabel('$B_n$','Interpreter','Latex');
set(gca,'FontSize',15);
subplot(1,2,2);
surf(KMesh,BMesh,squeeze((1./var_policy.Pf(z1,:,:))-1)*12);
title('Interest Rate');
xlabel('$K$','Interpreter','Latex');
ylabel('$B_n$','Interpreter','Latex');
set(gca,'FontSize',15);

As shown, given capital, as the bond holding of non-stock holders increases, the demand for bond increases and for stock decreases, which results in a lower interest rate, lower stock price and lower equity premium.
We can then simulate the economy
>> SimuRslt = simulate_guvenen2009(IterRslt);
Periods: 1000
shock K bn_shr Y c_h c_n Inv Ps Div Pfbn_shr_next Knext omega PDratio Rs R_epSharpeRatio
8 126.1 0.3264 4.267 5.103 3.003 0.8448 111.7 0.3939 0.9978 0.3265 126.1 0.9521 283.7 1.0040.002244 0.04976
Elapsed time is 28.379607 seconds.
Periods: 2000
shock K bn_shr Y c_h c_n Inv Ps Div Pfbn_shr_next Knext omega PDratio Rs R_epSharpeRatio
9 127.3 0.3652 4.341 5.197 3.04 0.869 119.3 0.4121 0.9989 0.3655 127.3 0.9493 289.5 1.003 0.00224 0.05007
Elapsed time is 25.144899 seconds.
Periods: 3000
shock K bn_shr Y c_h c_n Inv Ps Div Pfbn_shr_next Knext omega PDratio Rs R_epSharpeRatio
11 135.2 0.5826 4.547 5.381 3.158 0.9444 136.5 0.4421 1.001 0.5829 135.3 0.928 308.7 1.0010.002332 0.05232
Elapsed time is 14.973382 seconds.
and inspect the distribution of state variables in the ergodic set
startPeriods = 500;
Kss = IterRslt.params.Kss;
K = SimuRslt.K(:,startPeriods:end);
omega = SimuRslt.omega(:,startPeriods:end);
figure; clf;
subplot(1,2,1);
histogram(K/Kss,40,'Normalization','probability');%,'BinLimits',[0.5,0.75]);
title('Distribution of $K$','Interpreter','Latex');
xlabel('$K/K^{ss}$','Interpreter','Latex');
ylabel('Fractions');
set(gca,'FontSize',15);
subplot(1,2,2);
histogram(omega,40,'Normalization','probability');%,'BinLimits',[0.5,0.75]);
title('Wealth Share Distribution','Interpreter','Latex');
xlabel('Wealth of Stockholders','Interpreter','Latex');
ylabel('Fractions');
set(gca,'FontSize',15);

What’s Next?¶
This example illustrates how to work with parameterized expectation with the toolbox, which greatly reduces computation cost for models that can be solved with predetermined endogenous state variables. The finite-agent asset pricing model is an important class in the Macro and Macro-finance literature. Since the framework of the toolbox solves allocations and prices jointly respecting the optimality and market clearing conditions, it is able to deal with setups which are traditionally regarded as ill-conditioned, such as models with very high risk aversion levels (see example Barro et al. (2017)) or models with equilibrium asset returns that are close to collinear such as Cao, Evans, and Luo (2020).
Or you can directly proceed to Toolbox API.