Computing Equilibrium Bond Prices in the Vayanos-Vila Model

We develop tools for computing equilibrium bond prices for the discrete-time version of the Vayanos-Vila (2009) model. With the maturity structure included in pricing factors, factor loadings for equilibrium bond yields depends on parameters describing maturity structure dynamics and other model parameters. An illustrative example shows that the effect on the yield curve of a supply shock originating in a given maturity, although hump-shaped around the originating maturity, is to change yields broadly across all maturities.


Introduction and Summary
Governments often attempt to influence the yield curve through open-market operations. 1 This paper theoretically studies the effect of changes in the maturity structure engineered by such policy actions. It does so by developing tools for computing equilibrium bond prices for the explicit model of bond markets of Vayanos and Vila (2009).
The Vayanos-Vila model has two types of investors, risk-averse arbitrageurs and a set of maturity-specific preferred-habitat investors who differ from each other in the preferred maturity of bonds they hold. Equilibrium bond prices and yields are affine functions of pricing factors including supply factors affecting the net supply of bonds by preferred-habitat investors. If the model is specialized so that the net supply is price-inelastic, then the preferred habitat investors as a whole can be interpreted as the government exogenously supplying bonds of different maturities, and the maturity structure becomes a linear combination of the supply factors.
This specialized model has been put to use in recent empirical studies. The dynamics of the maturity structure in Greenwood and Vayanos (2014) and Greenwood et. al. (2015) is very restrictive: only one supply factor drives the whole profile of the maturity structure, even though (as remarked in Greenwood and Vayanos (2014, p. 12)) multiple factors might be needed to describe supply. D' Amico and King (2013), in their analysis of the recent Federal Reserve's large-scale purchase of Treasury securities, considered the supply effect of government bonds of a wide variety of maturities. To compute the factor loadings for equilibrium yields, they resort to an 1 Recent examples include the Federal Reserve's LSAP (Large-Scale Asset Purchases) and the Bank of Japan's QQE (quantitative and qualitative easing), both of which involve a very large amount of purchases of long-term government bonds. approximation suggested by Vayanos and Vila (2009). The paper by Hamilton and Wu (2012) is another empirical study of the supply effect motivated by the Vayanos-Vila model. However, their regressions linking bond supply to bond prices are not derived from the model. 2 The empirical applications are constrained thus because computing equilibrium bond prices in the Vayanos-Vila model with more-than-several pricing factors has been viewed as intractable. The difficulty lies in the system of equations for the affine factor loadings. Unlike in the standard ATSM (affine term structure model) in modern finance, it is a non-linear system and cannot be solved recursively.
The contribution of this paper is to provide tools for overcoming this computational difficulty, for the discrete-time version of the Vayanos-Vila model. The key observation is that the nonlinear factor loadings equations are QVEs (quadratic vector equations). We show that the solution can be obtained via a system of ODEs (ordinary differential equations) in a single variable. We verify that, for a set of model parameters we examine, the solution thus obtained can be located quickly by a fixed-point algorithm, a finding which should prove useful in a structural estimation of the model.
Including a large number of pricing factors in the affine yield equations, made feasible by the tools just described, is consistent with the well-known empirical fact that only up to three factors are enough to describe the yield curve. This is because the affine factor loadings, being a 2 There are two reasons, one technical and the other substantive. The technical reason is that they ignore a connection implied by the model between factor loadings for yields and those for risk premia. The substantive reason is that they do not include bond supply in the factors affecting bond prices. Their inclusion of bond supply in their regressions explaining bond prices is therefore not justified by the model. Cochrane (2011) puts this latter point more bluntly: "Hamilton and Wu do not use the term structure model to infer the effect of bond supply. The results are just a simulation from a regression of yield changes on bond supply." (Italics original) solution to the model, are not free parameters. The model's parameters are the VAR (vector autoregression) parameters describing maturity structure dynamics and a risk-aversion parameter for the arbitrageurs. Therefore, given the VAR parameters, the factor loadings depend on only one parameter. In this sense, the model is a one-factor model.
As an illustration, we compute equilibrium bond prices for the case of maximal number of supply factors. That is, the maturity structure itself is the vector of supply factors, or there are as many supply factors as there are maturities. This is an interesting case to study, for two reasons.
First, the maturity structure dynamics can account for what we call the "legacy" property of bonds, which is that an increase in the supply of n-period bonds necessarily leads to an increase, in the next period, of the same amount in the supply of (n − 1)-period bonds, unless there is an offsetting market operation. Second, we can give a rigorous formulation of a "local" supply shock. Shocks to the maturity structure, one for each maturity, are local in that they are uncorrelated across maturities.
Our computation shows that he effect of a positive local supply shock originating at a given maturity is to shift the yield curve up. Although hump-shaped around the originating maturity, the shift occurs broadly across all maturities. This finding is surprising given the conjecture by Vayanos and Vila (2009) that the effect of a local shock would be local in the yield curve. 3 The rest of the paper is organized as follows. Section 2, in order to set a stage for our computational tools in Section 3, presents the specialized Vayanos-Vila model in discrete time.
Section 4 examines the effect of local supply shocks for the illustrative example. Section 5 states the conclusion and an agenda for future research.

The Specialized Vayanos-Vila Model in Discrete Time
The bond pricing model to be considered in this paper is the discrete-time version of Vayanos and Vila (2009) with the added restriction that the preferred-habitat investors are replaced by an entity interpretable as the government supplying bonds inelastically. This section presents the model and a system of equations for equilibrium bond prices. The equation system will be the object of analysis in the next section.

The Discrete-Time Bond Pricing Model
Denote by P (n) t the price of zero-coupon government bonds with maturity n in period t. The (continuously compounded) yield to maturity y (n) t is related to the bond price P (n) t as in (2.1) The nominal share of the n-period bond supplied is denoted by s (n) t . The longest maturity of bonds supplied is N. The shares add up to one: N n=1 s (n) is the maturity structure of bonds supplied.
Those bonds are traded by arbitrageurs. If z (n) t is the nominal share of their n-period bond holdings, the holding-period return on their bond portfolio is The decision problem of arbitrageurs is to maximize the risk-adjusted portfolio return subject to the adding-up constraint: 4 Here, γ (≥ 0) is a risk aversion coefficient. The FOCs (first-order conditions), derived in Appendix 1 for the sake of completeness, are which state that, for arbitrageurs to hold a given amount of n-period bonds, the risk premium on those bonds must be enough to compensate for an increase in the volatility of the portfolio return R t+1 due to the inclusion of those n-period bonds in the portfolio.
The bond market equilibrium is that z (n) t = s (n) t for n = 1, 2, ..., N. Because of the adding-up constraint N n=1 z (n) t = 1, only N − 1 of these equilibrium conditions are independent.
The choice (to be mentioned right below) of factors for bond pricing implies that the price of one-period bonds, P (1) t , is determined outside the model. Given P (1) t , the N − 1 equilibrium The equilibrium bond prices are assumed to be a time-invariant function of a vector, f t , of factors whose dynamics is given by a Gaussian VAR(1): (2.5) The shortest yield y (1) t is included in f t , which means that the short rate, and hence the price of one-period bonds P (1) t , are given to the model. The remaining factors are denoted by the (2.6) We will call the factors included in β t the supply factors because they span the maturity structure . (2.8)

The Factor Loadings Equations
Vayanos and Vila (2009) have shown for the continuous-time case that the time-invariant function for equilibrium bond prices is exponential affine, that is, log P (n) t =ā n +b n f t , so by (2.1), y (n) t = a n + b n f t where a n ≡ −ā n n , b n ≡ −b n n , n = 1, 2, ..., N. (2.9) They also derived a system of equations that have to be satisfied by the affine factor loadings (i.e., the coefficients (ā n ,b n ) in the above affine function).
For the discrete-time model at hand, the affine functional form remains valid. The factor loadings equations, derived in Appendix 1 for completeness's sake, arē where F is the number of factors (which equals 1 + K where K is the number of supply factors).

Yields and Risk Premia
For later use in Section 4, we point to a formula showing how yields are connected to risk premia.
The usual forward iteration argument (included in Appendix 2) then yields: where A n ≡ − 1 2 n−2 k=0b n−k−1 Ωb n−k−1 is the Jensen inequality term. This formula shows that the yield on an n-period bond has two time-varying components. One is the expectations hypothesis benchmark, namely the average over the bond's life of the current and future short rates. The other component, which is sometimes called the term premium, can be written as the average of the longitudinal profile of the risk premium, rp (n) t , rp (n−1) t+1 , ..., rp (2) t+n−2 , on the n-period bond as it approaches maturity.

Solving the Quadratic Vector Equation (2.10)
The QVE Unlike in the standard textbook ATSM (affine term structure model) and in contrast to (2.11), equation (2.10) is not a recursion thanks to the quadratic term (the second term on the right-hand side of the equation) that involves the whole array (b 1 , ...,b N−1 ). To make it clear that it is a QVE (quadratic vector equation), define the NF-dimensional stacked vectorb as (recall: F is the (3.1) Then it is straightforward to show that (2.10) and the initial conditionb 1 = δ ≡ (−1, 0, ..., 0) together can be written as 5 Incidentally, there is a corresponding QVE for risk premia. Appendix 1 shows that risk premia are linear in the factor vector (i.e., affine with zero intercepts). If h n is the vector of factor loadings for the risk premium on bonds of maturity n, It is related to the bond price factor loadingsb n as (see Form the NF × 1 stacked vector h of (h 1 , h 2 , ..., h N ). By (5), h is related tob as Solving this forb and substituting into (3.2), we obtain the QVE for the risk-premium factor loadings: . (3.6)

The ODE
If arbitrageurs are risk neutral, i.e., if the risk aversion parameter γ in (2.3) is zero, then the equation (2.10) or (3.2) has a unique solution Although the square matrix M, which is of size NF, can be large, computingb * is fast because there is a well-known analytical expression for the inverse of the large matrix M. 6 If γ > 0, because of the quadratic term g(b), there can be multiple solutions if a solution exists. We seek 6 Here is the expression for M −1 : the solution that converges tob * as γ ↓ 0. 7 To study such a solution for a given value γ > 0 of the risk aversion parameter, write the equation (3.2) as (3.9) By construction,b =b * is the only solution to (3.9) whenγ = 0. Assume that the model parameters are such that there is one or more solutions over an intervalγ ∈ (0, U). 8 Since f(., .) is continuously differentiable in bothb andγ, each solution is a differentiable function ofγ by the Implicit Function Theorem. For any given solution, write the function asb(.): (by the definition in (3.9) of f(b, γ)). (3.10) The functionb(.) that satisfies the conditionb(γ) →b * asγ ↓ 0 can be obtained by viewing (3.10) as an ODE (ordinary differential equation) inγ and solve it on the interval [0, γ] with the 7 As discussed in McCafferty and Driskill (1980), multiple equilibria is generic to nonlinear rational expectations models. To cite Greenwood and Vayanos (2014): "If yields are highly sensitive to shocks to the supply risk factor..., then bonds become highly risky for arbitrageurs. Hence, arbitrageurs absorb supply shocks only if they are compensated by large changes in yields, making the high sensitivity of yields to shocks self-fulfilling." They rule out such equilibria and note that the equilibrium they select "is well behaved in the sense that when a (their risk aversion parameter) converges to zero, it converges to the unique equilibrium that exists for a = 0." initial conditionb(0) =b * . The value of the solution path atγ = γ is the desired solution to the QVE for a given value γ of the risk aversion parameter.

An Iterative Fixed-Point Algorithm
Albeit not requiring a numerical calculation of the Jacobian matrix ∂g(b(γ)) ∂b , this method can be computationally demanding if the dimension of the matrix to be inverted, M +γ An obvious alternative is to utilize the fixed-point iteration starting fromb (0) =b * . This alternative method, with no need to do numerical matrix inversion, is faster by orders of magnitude. We have not been able to characterize a condition under which this iteration converges, though. 9 However, for numerical examples we experimented and certainly for the example of the next section, this iteration quickly converges to the solution obtained via the ODE (3.10).

An Illustrative Example
The model parameters are: γ (the arbitrageurs' risk aversion parameter in (2.3)), N (the maximum maturity), (c, Φ, Ω) (the VAR parameters describing the factor dynamics for the factor vector f t , see (2.5)), and C (the matrix connecting supply factors β t to the maturity structure (s (2) t , ..., s (N) t ), see (2.7)). The factor vector f t consists of the short rate and the supply factors β t , as in (2.6). The previous section was about how to compute the affine factor loadings (the coefficients of f t , (ā n ,b n ), n = 1, 2, ..., N, in (2.9)) given the model parameters.
With computation tools for an arbitrary number of supply factors in hand, we can allow the matrix C in (2.7) to be of full rank. Without loss of generality, we can set C = I N−1 . So the supply factors β t are the maturity structure itself: (4.1) Consequently, the supply shocks (ε (2) t , ε (3) t , ..., ε (N) t ) are now shocks to the maturity structure, K (the number of supply factors) equals N − 1, and F (the dimension of the factor vector f t ) equals 9 The sufficient condition identified in Poloni (2013, Theorem 4) requires that both Φ and Ω be nonnegative.
N (= 1 + K). The affine function for yields in (2.9) is now If b jn is the j-th element of b n , the profile (b j1 , b j2 , ..., b jN ) represents the impact response of the entire yield curve to a one-unit shock originating in maturity j = 2, 3, ..., N. That is, In this section, we consider a specific example of the factor dynamics with this choice (4.1) of the factor vector. The example's basic features are: (i) supply shocks are purely local in the maturity spectrum in that they are uncorrelated across maturities, (ii) shock variances are the same for all maturities (so maturity-specific effects on yields, if any, are not attributable to a supply volatility unique to the maturity), (iii) time dependence in the maturity structure is due solely to what we call the "legacy" nature of bonds that an n-period bond becomes an (n − 1)-period bond in the next period, and (iv) the short rate is exogenous to the maturity structure (so the supply effect on the yield curve is not contaminated by the expectations effect that operates through changes in the future path of the short rate).

Specifying the Factor Dynamics
The VAR (2.5) is specified as 10 (short rate) y (1) where the bond supply shocks (ε 2t ,...,ε Nt ) are uncorrelated across maturities. The short rate y (1) t is exogenous to the maturity structure, because the short rate shock ε 1t is uncorrelated with the bond supply shocks and also because the lagged supply of bonds do not enter the short rate equation.
The extent of the bond "legacy" property is represented by θ, the fraction of the n-period bond supply that will be left as (n − 1)-period bonds in the next period. It equals 0 if there is a fully offsetting market operation in the next period, and 1 if there are no offsetting operations. 10 The implied restrictions on the VAR parameters (Φ, Ω) are as follows.

Calibration
To avoid heavy demand on CPU time, we choose the unit interval to be a quarter, not a month, and set N to 80 quarters (20 years). 11 So y (1) t is the 3-month interest rate. For calibration purposes we use the U.S. Treasury zero-coupon yield curve calculated by Gurkaynak, Sack, and Wright (2007). 12 The factor dynamics consists of the AR (1)  Regarding the maturity structure dynamics (4.5), we pick the intercept terms (c 2 , ..., c N ) so t for any maturity n. We take the yield for the quarter to be for the last business day of the quarter. For the short rate we use the value implied by the parametrized yield curve at n = 1 (3 months). This short rate series is very similar to the constant-maturity 3-month rate available from Federal Reserve Board's Table H -15. that the steady-state value of s (n) t is 1/N for n = 2, ..., N. 13 We set σ n = λ/N with λ = 0.01 for n = 2, 3, ..., N, so the volatility of the supply shock is uniform across maturities. The results to be reported below are not sensitive to the choice of λ over the range [0.001, 0.1].
This leaves two parameters, γ (the risk aversion parameter) and θ (the legacy parameter that controls the maturity structure dynamics). The latter has no effect on the steady-state yield curve. 14 We can therefore calibrate γ to the steady-state yield curve without specifying the value of θ. We set γ = 48 because, as shown in Figure 1,

Response of Risk Premia to Local Supply Shocks
We are now ready to describe the effect of bond supply shocks. We first examine the effect on risk premia before turning to yields, in order to exploit the formula (2.13) connecting risk premia to yields. This is also the route taken by Greenwood and Vayanos (2014) for their model of a single supply factor. It is straightforward to show (see (A1.3)) that risk premia are linear (affine with 13 Lets t t ) be the maturity structure,c ≡ (c 2 , ..., c N−1 , c N ) , and letΦ be as in footnote 10. Then the steady-state value ofs t , call its, satisfies: (I −Φ)s =c. The factor loadings equations (2.10) for {b n } do not involve the VAR intercepts c ≡ (c 1 , c 2 , ..., c N ). Therefore, the maturity structure intercepts (c 2 , ..., c N ) do not affect the response profiles of yields and risk premia to supply shocks, which are a function of (b 1 , ...,b N ), to be displayed in Figures 2 and 3. They do affect the intercepts (a 1 , a 2 , ..., a N ) in the yield affine function in (2.9) and also the factor vector f t .
zero intercepts) in f t . Let h jn be the coefficient of the j-th element of the coefficient vector h n . 16 The impact response of the term structure of risk premia to a local supply shock originating in maturity j ( j = 2, 3, ..., N) is given by the profile (h j1 , h j2 , ..., h jN ). As it turns out, the factor loadings for risk premia, and hence the impact response profile, are invariant to θ. 17 This θ-invariant impact response profile for risk premia is shown in Figure 2 for three originating maturities, 20 quarters (5 years), 40 quarters (10 years), and 80 quarters (20 years).
The shock is a 1 percentage point increase in the share of bonds of the originating maturity. There are two notable features: • The profile shifts up as the originating maturity gets longer. That is, ∂ε jt is increasing in j for each n = 2, 3, ..., N. An increased supply of bonds of maturity j (= 2, 3, ..., N), which is to be purchased by arbitrageurs, makes the bond portfolio held by arbitrageurs more sensitive to changes in the short rate. Thus the portfolio's volatility (the Var(R t+1 ) in the first-order conditions (2.4)) increases, which raises the risk premium on bonds of any given maturity. This effect is stronger if the originating maturity (j) of the additional bonds supplied is longer, because the portfolio's volatility increase is greater.
• The profiles are upward-sloping. That is, ∂ε jt is increasing in n as well. The risk premium needed to compensate arbitrageurs for the portfolio's volatility increase is greater for longer-term bonds in the portfolio because they are more sensitive to changes in the short rate. 18 16 Footnote 5 shows how {h n } are related to the yield factor loadings {b n }.
17 This property holds because the only source of time dependence in the maturity structure is the legacy property. It does not hold if the diagonal elements of the VAR coefficient matrixΦ shown in footnote 10 are not zero.

Response of the Yield Curve to Local Supply Shocks
Differentiating both sides of the formula (2.13) with respect to ε jt and noting that (i) the shock does not affect current and future short rates and (ii) subsequent responses, ∂ rp (n) t+k ∂ε jt (k = 1, 2, ...), are known at date t, we obtain the formula that links risk-premium responses to the yield response: (4.7) That is, the impact response of the n-period yield to a supply shock originating in maturity j is the average of the longitudinal profile of risk-premium responses over the life of the n-period bond.
The longitudinal risk-premium profile is related to the cross-section risk-premium profile of impact responses, according to the following formula: is increasing in j.
(b) The longitudinal profile is a declining sequence, because θ ≤ 1 and because ∂ rp (n) t ∂ε jt declines Greenwood and Vayanos (2014, Proposition 2). with both j and n (i.e., ∂ rp (n) t ∂ε jt is increasing in both j and n, as noted above). • For each originating maturity, the profile is greater for risk premia than for yields. That is, ∂ε jt for all n, j = 2, 3, ..., N. Even with θ = 1, the profiles shown in Figure 2 lie above the corresponding profiles shown in Figure 3 for each originating maturity. The reason is that , whose average equals ∂y (n) t ∂ε jt by (4.7), is a declining sequence, as noted in (b) above. 19 • As in Figure 2 for risk premia, the yield profile shifts up as the originating maturity gets longer.
• In contrast to the profiles shown in Figure 2, the yield profiles are hump-shaped around the originating maturity. The reason is that the effect of a supply shock at originating maturity j lasts only up to j − 1 periods. For example, ε 2t ceases to exert yield-curve effects in the next period t + 1 when two-period bonds turn into one-period securities. For bonds whose maturity n is greater than j, the longitudinal profile in (4.8) is truncated, as in (4.9) There are two opposing effects on the average of this profile as n gets longer. One is the dilution effect due to the truncation after j − 1 periods. The other is that the profile before truncation is higher because the risk-premium effect increases with n. The former effect dominates the latter as n gets past j. Appendix 1. Derivation of Factor-Loading Equations (2.10) and (2.11)

Conclusion and Agenda for Future Research
We proceed in two steps: (a) Under the VAR dynamics (2.5) and the conjecture of the affine bond pricing (2.9), derive the arbitrageurs' FOCs (first-order conditions) that involve the affine coefficients (ā n ,b n ).
(b) Impose on the FOCs the bond market equilibrium conditions, z Regarding (a), the FOCs for n-period bonds, reproduced from the text and to be derived below, are We can rewrite the risk premium (the left-hand side of (A1.1)) and 1 2 ∂ Var t (R t+1 ) ∂z (n) t (on the right-hand side of (A1.1)) as Turning to (b), upon the imposition of the market equilibrium conditions, the expression (A1.4) for can be written as (by (2.7)) (A1.5) Substituting (A1.3) and (A1.5) into (A1.1), we obtain a set of equations that is affine in f t . This has to hold for any f t . Setting both the constant term and the coefficients of f t to zero, we obtain the system of equations, (2.10) and (2.11) in the text, for affine factor loadings (ā n ,b n ).
• Show: for n = 2, 3, ..., N. This can be obtained easily from two well-known formulas. One is that, for a normally distributed random variable X, E[exp(X)] = exp[E(X) + 1 2 Var(X)], and the other is exp(x) ≈ 1 + x for x ≈ 0. In the former, we can set X = p (n−1) t+1 − p (n) t because, as shown above, p (n−1) t is normally distributed conditional on date t information. In the latter formula, set • Show: Var t (R t+1 ) ≈ d t Ωd t where d t ≡ Note: The impulse response of the term structure of risk premia. The shock is a 1 percentage point increase in the share of bonds of the originating maturity.