Symbolic stationarization of dynamic equilibrium models

Canova, Fabio; Paulsen, Kenneth Sæterhagen Working Paper Symbolic stationarization of dynamic equilibrium models Working Paper, No. 18/2021 Provided in Cooperation with: Norges Bank, Oslo Suggested Citation: Canova, Fabio; Paulsen, Kenneth Sæterhagen (2021) : Symbolic stationarization of dynamic equilibrium models, Working Paper, No. 18/2021, ISBN 978-82-8379-217-1, Norges Bank, Oslo, https://hdl.handle.net/11250/2835495 This Version is available at: https://hdl.handle.net/10419/264937 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. https://creativecommons.org/licenses/by-nc-nd/4.0/deed.no Symbolic Stationarization of Symbolic Stationarization of Dynamic Equilibrium ModelsDynamic Equilibrium Models NORGES BANK RESEARCH 18 | 2021 Fabio Canova Kenneth Sæterhagen Paulsen WORKING PAPER NORGES BANK WORKING PAPER XX | 2014 RAPPORTNAVN 2 Working papers fra Norges Bank, fra 1992/1 til 2009/2 kan bestilles over e-post: [email protected] Fra 1999 og senere er publikasjonene tilgjengelige på www.norges-bank.no Working papers inneholder forskningsarbeider og utredninger som vanligvis ikke har fått sin endelige form. Hensikten er blant annet at forfatteren kan motta kommentarer fra kolleger og andre interesserte. Synspunkter og konklusjoner i arbeidene står for forfatternes regning. Working papers from Norges Bank, from 1992/1 to 2009/2 can be ordered by e-mail: [email protected] Working papers from 1999 onwards are available on www.norges-bank.no Norges Bank’s working papers present research projects and reports (not usually in their final form) and are intended inter alia to enable the author to benefit from the comments of colleagues and other interested parties. Views and conclusions expressed in working papers are the responsibility of the authors alone. ISSN 1502-8190 (online) 978-82-8379-217-1 (online) Symbolic Stationarization of Dynamic Equilibrium Models Fabio Canova * Kenneth Sæterhagen Paulsen  December 20, 2021 Abstract Dynamic equilibrium models are specied to track time series with unit root-like behavior. Thus, unit roots are typically introduced and the optimality conditions adjusted. This step requires tedious algebra and often leads to algebraic mistakes, especially in models with several unit roots. We propose a symbolic algorithm that simplies the step of rendering non-stationary models stationary. It is easy to implement and works when trends are stochastic or deterministic, exogenous or endogenous. Three examples illustrate the mechanics and the properties of the approach. A comparison with existing methods is provided. Keywords: DSGE models, unit roots, endogenous growth, symbolic computation. * BI Norwegian Business School and CEPR, [email protected]  Norges Bank, corresponding author; [email protected]. This working paper should not be reported as representing the views of Norges Bank. The views expressed in this paper are those of the authors and do not necessarily reflect those of Norges Bank. 1 1 Introduction Dynamic Stochastic General Equilibrium (DSGE) models are popular in academic and policy institutions and used for a variety of purposes: informally matching stylized facts present in the data, structural estimation, scenario analyses, and forecasting. Often these models just consider stationary disturbances and a standard production structure with constant returns to scale, making the output a vector of stationary time series. When this is the case, one has to choose how to compare the model to data displaying near non-stationary behavior and how to design estimation and inferential procedures which can cope with the mismatch (see e.g. Canova (2014)). Alternatively, one could assume that some of the disturbances are non-stationary, for example, assume permanent shocks to total factor productivity. In this case, a balanced growth path must be found prior to the computation of the solution and the optimality conditions transformed so that a stationary equilibrium exits. Standard solution procedures such as Blanchard and Kahn (1980), Klein (2000) or Sims (2010), in fact, work under the assumption that the variables entering the optimality conditions are stationary. Manually performing the stationarization step is cumbersome, it involves tedious algebra and, more importantly, it is prone to mistakes, especially in medium or large scale models featuring a number of non-stationary disturbances. This paper presents a fast and easy-to-implement algorithm that mechanically performs the stationarization step in models featuring unit roots or deterministic trends, when the sources of the trends are exogenous or endogenous, such as in Comin and Gertler (2006). Because the approach is symbolic, the cost of computations is trivial. Thus a researcher may explore the consequences of alternative assumptions regarding the location of the unit roots (preferences, technologies, policy shocks), examine the equations that are aected, and the implications they have for the observable variables. Such an exercise allows one to improve the specication and the empirical content by confronting the long run implications of the model and of the data systematically. The output of the algorithm consists of a system of equations containing the stationary version of the original non-stationary model, and it can be written to a text le or exported to LaTeX. The user can thus check that the solution makes sense and then perform a standard dynamic analysis. The algorithm is contained in a package of routines, called NB toolbox (Paulsen (2021)), which also has a number of ready-built functions a researcher can employ for standard dynamic analysis. Given that the stationary model can be written to a text le, the output of the algorithm can also be used as an input in other toolboxes, such as DYNARE, Dynare (2021); IRIS, GPMN (2020); or RISE, Maih (2021), for ltering or estimation, if the user wishes to do so. IRIS and DYNARE have the capability to deal with non-stationary disturbances in dynamic models. Our approach is more general than both. The rst toolbox uses a numerical approach to perform the 2 stationarization step. Thus, the solution features approximation errors. In addition, the user needs to provide the balanced growth path for the approximation to be performed. The second toolbox needs the user to provide both the growth factors and the variables which are associated with the growth factors. In addition, at the moment, it can handle only exogenous unit root processes. The rest of the paper is organized as follows. The next section describes the problem to be solved. The algorithm is in Section 3. Section 4presents the syntax of the model le and the commands. Section 5 has examples. In Section 5.1 we consider a RBC model driven by non-stationary labor augmenting technology shocks. We stationarize the model manually (Section 5.1.1); with our algorithm (Section 5.1.2) and compare the two solutions. In Section 5.2 we repeat the exercise for a larger-scale model with unit root in the labor-augmenting and the investment-specic technologies, as in Justiniano et al. (2011)). We compare its long run properties to those of a model where unit roots appear in the production technology and in preferences. Section 5.3 considers a model with endogenous growth and compares manual and automatic solutions. Section 6discusses the relationship between our approach and existing methods. Section 7concludes. A number of on-line appendices contain the details and the les mentioned in the text. 2 The problem The optimality conditions of a DSGE model can be written as: Et[F(Xt+1, Xt, Xt−1, At, At−1, Ut, θ)] = 0, (1) where F(.) is a M×1 vector of continuous non-linear functions, and Xt is a M×1 vector of endogenous variables, with the convention that the rst R variables are non-stationary, and the remaining M-R are stationary, and Ut is a N×1 vector of stationary exogenous forcing variables. At is a Q×1 vector of non-stationary forcing variables which, for the sake of the presentation, are assumed to enter without leads, even though this is not a limitation of the approach, as leads can always be considered by extending the vector At . At may be exogenous or endogenous and, in this latter case, may depend on a subset of Ut, θ and Xt . Finally, θ is a p×1 vector of structural parameters. The presence of the At in (1) rules out the possibility of using standard perturbation methods to solve the problem, since the Blanchard and Kahn condition fails. However, if the problem has a stationary representation, we can redene the R endogenous non-stationary variables as Zi t=Xi t Hi(At),∀i∈[1, R] and the M-R stationary variables as Zi t=Xi t,∀i∈[R+ 1, M] , and proceed in the standard way. Here Hi are continuous dierentiable functions of At and represents the balanced growth path of the model. Let ∆Zi t≡Hi(At) Hi(At−1) Xi t−1 Xi t ∀i∈[1, R] and ∆At=At At−1 and set Xt−1 Xt= 1 . Rather than nding the Hi(At) functions, which is complicated when R is large, it is easier to compute the functions Gi satisfying Gi(log(∆Zt), log(∆At)) = 0,∀i∈[1, R] . These functions also characterize the 3 balanced growth path and are obtained with a set of restrictions, stated in Section 3.1. If we let Yt=  Zt ∆Zt ∆At , (2) a stationary representation of the optimality conditions of the model is Et[ˆ F(Yt+1, Yt, Yt−1, Ut, θ)] = 0, (3) where ˆ F(.) is a function representing the non-linear equations, and Yt has dimension M+R+Q . Note that ∆Zt and ∆At belong to a block exogenous system when At is exogenous. In words, our algorithm creates the Yt from any set of (Xt, At, Ut) , and transforms (1) into (3). To perform the transformation one needs four steps: 1. Find the restrictions holding in the balanced growth path. 2. Identify the R non-stationary variables. 3. Compute the functions Gi,∀i∈[1, R] . 4. Rewrite (1) into (3). Clearly, trends must be properly placed for the procedure to succeed, as the restrictions obtained from the rst step need not be consistent with balanced growth. For example, it is well know that, with a unit root in government expenditure, there are no restrictions on preferences and technologies producing balanced growth in a standard real business cycle (RBC) model. 3 The algorithm This section describes the algorithm that automatizes the four steps outlined in Section 2. 3.1 Finding the restrictions holding in the balanced growth One method to quickly identify which variables are non-stationary is to nd the set of restrictions that must hold to rewrite the model in a stationary form. These restrictions are generically given by: K(X∆ ss, A∆ ss, θ) = 0, (4) where X∆ t=log(Xt Xt−1) and A∆ t=log(At At−1) , and X∆ ss and A∆ ss are the values along the balanced growth path. Here the function K represents a vector of L linear restrictions. In general, L≥M . Let θp and θd be two parameters, h and g two generic functions and h± a function formed using h±g . Let c and 4 v be functions obtained applying the rules described below to h and g , respectively. We can nd the restrictions given by K , by applying the following rules to the equations of (1) by forward accumulation, and eliminating all duplicated restrictions 1 . Equal (=)  If Xss =Yss , then X∆ ss =Y∆ ss .  If Xss =h(Yss, Zss) , then X∆ ss =c(Y∆ ss , Z∆ ss) .  If Xss =h±(Yss, Zss) , then X∆ ss =c(Y∆ ss , Z∆ ss) and X∆ ss =v(Y∆ ss , Z∆ ss) .  If h(Yss, Zss) = θp then c(Y∆ ss , Z∆ ss)=0 . To understand the meaning of these restrictions, note that the rst one states that if a model has an equation of the form Xss =Yss , Xss and Yss need to grow at the same rate along the balanced growth path. Similarly, the second restriction forces the two sides of the equation to have the same growth rate. The third expression consists of two terms; Xss =h(Yss, Zss)±g(Yss, Zss) and forces them to grow at the same rate along the balanced growth path. The last condition requires h(Yss, Zss) to display no growth, as θp is a constant. The logic of the next expressions are similar: Plus/Minus ( ± )  If Xss ±Yss , then X∆ ss =Y∆ ss .  If Xss ±Yss ±Zss , then X∆ ss =Y∆ ss and X∆ ss =Z∆ ss .  If Xss ±h(Yss, Zss) , then X∆ ss =c(Y∆ ss , Z∆ ss) .  If Xss ±h±(Yss, Zss) , then X∆ ss =c(Y∆ ss , Z∆ ss) and X∆ ss =v(Y∆ ss , Z∆ ss) .  If h(Yss, Zss)±θp , then g(Y∆ ss , Z∆ ss) = 0 .  If θp±θd , then θp±θd . The last rule allows us to keep track of parameters. Since a parameter is an object of an equation, we need to ensure that, when applying the relevant operators, they are properly accounted for. Times (*)  If XssYss , then X∆ ss +Y∆ ss .  If XssYssZss , then X∆ ss +Y∆ ss +Z∆ ss .  If Xssh(Yss, Zss) , then Xss +c(Y∆ ss , Z∆ ss) . 1 In a forward accumulation exercise, one rst xes the independent variable with respect to which the scaling is performed and then computes the trend of each sub-expression recursively. 5  If Xssh±(Yss, Zss) , then c(Y∆ ss , Z∆ ss) = v(Y∆ ss , Z∆ ss) and X∆ ss +c(Y∆ ss , Z∆ ss) .  If Xssθp , then X∆ ss .  If θp·θd then θp·θd . The last two rules imply that multiplying Xss by a parameter leaves its growth rate unchanged, and that parameters do not grow. Divide (/)  If Xss Yss , then X∆ ss −Y∆ ss .  If (Xss Yss ) Zss , then X∆ ss −Y∆ ss −Z∆ ss .  If Xss (Yss Zss ) , then X∆ ss −Y∆ ss +Z∆ ss .  If Xss h±(Yss,Zss) , then c(Y∆ ss , Z∆ ss) = v(Y∆ ss , Z∆ ss) and X∆ ss −c(Y∆ ss , Z∆ ss) .  If Xss θp , then X∆ ss .  If θp θd then θp θd . Power ()  If (Xss)θp then θp∗X∆ ss .  If θYss p , then Y∆ ss = 0 , i.e. Yt must be stationary.  If θh(Yss,Zss) p , then c(Y∆ ss , Z∆ ss) = 0 , i.e. h(Yt, Zt) must be a stationary.  If θθd p then θθd p . Expression of the form (Xss)Yss or (Xss)h(Yss,Zss) are not supported. This is because, for example, the rst expression leads to a combination of two restrictions Y∆ ss = 0 and YssX∆ ss that prevent us from solving for the balanced growth path separately from the stationary steady state, as Yss appears in the system in (4). Situations of this type however are rare in the current wave of DSGEs. Exponential (exp)  If eXss then X∆ ss = 0 , i.e. Xt must be stationary.  If eh(Yss,Zss) then c(Y∆ ss , Z∆ ss)=0 , i.e. h(Yss, Zss) must be a stationary.  If eθp then eθp . Natural logarithm (log)  If log(Xss) then X∆ ss = 0 , i.e. Xt must be stationary. 6 reporting A = exp(cumsum(log(dA) - log(steady_state(dA)))); % level variables C = c*A; I = i*A; K = k*A; Y = y*A; c_log_dev = exp(log(c) - log(steady_state(c))); % deviations from steady state i_log_dev = exp(log(i) - log(steady_state(i))); k_log_dev = exp(log(k) - log(steady_state(k))); r_dev = r - steady_state(r); y_log_dev = exp(log(y) - log(steady_state(y))); With this le (saved as non_stationary.nb), the NB commands needed to stationarize the system and to compute impulse responses: modelNS = nb_dsge('nb_file','non_stationary.nb'); %Read the model file modelNS = set(modelNS,'name','Stationarize automatically'); % Name it param = struct( ); % Assign parameter values param.g = 1.03; param.gamma = 0.60; param.delta = 0.10; param.beta = 0.97; param.lambda = 0; param.std_u = 0.01; modelNS = assignParameters(modelNS,param); % Find the balanced growth path, (steps 1−3 of the algorithm): modelNS = solveBalancedGrowthPath(modelNS); % Make the model stationary (step 4 of the algorithm): modelNS = stationarize(modelNS) % Find the steady state and return it as a cell matrix: modelNS = checkSteadyState(modelNS,... 'solver','fsolve',... 'steady_state_solve', true,... 'steady_state_default', @ones); ss = getSteadyState(modelNS) % Return the balanced growth path solution as a cell matrix: bgp = getBalancedGrowthPath(modelNS) % Obtain the linear solution: y(t) = A y(t−1) + B eps(t) % and compute impulse responses modelNS = solve(modelNS); % Compute and plot impulse responses % Note modelS object represents the model stationarized manually. [irfs,~,plotter] = irf([modelNS,modelS],... 'variables',[modelNS.dependent.name(1:7),modelNS.reporting(:,1)'],... 'settings',{'legBox','off','legFontSize',18,'subPlotSize',[4,3],... 'figureTitle',false}); nb_graphInfoStructGUI(plotter) 13 The stationary model le produced by these commands is in Appendix B. If a user wants to save the stationary model for later purposes or to use in other packages, it is possible to write it to a le using: writeModel2File(modelNS,'stationarized.nb') We plot the dynamics induced by a one standard deviation increase in the labor-augmenting shock in the manual and automated solutions in Figure 1; the top part has the stationary variables, the bottom part the level variables. If the automated solution is accurate, the responses should be identical. y 5 10 15 20 25 30 35 40 1.82 1.822 1.824 1.826 1.828 1.82 1.822 1.824 1.826 1.828 r 5 10 15 20 25 30 35 40 0.0618 0.062 0.0622 0.0624 0.0626 0.0618 0.062 0.0622 0.0624 0.0626 l 5 10 15 20 25 30 35 40 0 1 2 0 1 2 k 5 10 15 20 25 30 35 40 4.61 4.62 4.63 4.64 4.65 4.61 4.62 4.63 4.64 4.65 i 5 10 15 20 25 30 35 40 0.5873 0.5874 0.5875 0.5873 0.5874 0.5875 dA 5 10 15 20 25 30 35 40 1.03 1.035 1.04 1.03 1.035 1.04 c 5 10 15 20 25 30 35 40 1.234 1.236 1.238 1.24 1.234 1.236 1.238 1.24 A 5 10 15 20 25 30 35 40 1 1.005 1.01 1 1.005 1.01 C 5 10 15 20 25 30 35 40 1.24 1.245 1.25 1.24 1.245 1.25 I 5 10 15 20 25 30 35 40 0.588 0.59 0.592 0.594 0.588 0.59 0.592 0.594 K 5 10 15 20 25 30 35 40 4.66 4.68 4.7 4.66 4.68 4.7 Y 5 10 15 20 25 30 35 40 1.83 1.84 1.85 1.83 1.84 1.85 AUTOMATIC STATIONARY MANUAL STATIONARY Figure 1: Impulse responses. A is the labour augmenting process. y, r, l, k, i, dA and c are the stationary variables in log deviation from the steady state; C, I, K and Y are the non-stationary levels of the variables. The two impulse responses are clearly indistinguishable. As the economy becomes permanently more productive, output (Y), consumption (C), investment (I) and capital (K) jump on impact and slowly reach the new long run level. Because the expected output to capital ratio increases on impact, the real interest rate (r) rises and this induces a fall in transitory consumption (c). As the economy moves to the 14 new steady state the real interest rate falls and transitory variables return to their original levels. 5.2 A New Keynesian model This section compares the manual and the automated decision rules of a slightly modied version of the Justiniano et al. (2011) model and studies the consequences of assuming unit roots in dierent shocks. The model allows us to illustrate how the algorithm works when the trend path is driven by multiple disturbances. The model deviates from Justiniano et al. (2011) in three ways: we assume Rotemberg rather than Calvo prices and wages adjustments; we use dierent functional forms for capital utilization and for investment-adjustment costs; and we set government purchases of goods and services to zero. There are six types of agents: nal and intermediate good producing rms, consumers, investment and capital producers, and entrepreneurs. Their problems and the optimality conditions are in the Appendix C. The model features nine disturbances, two of which are assumed to be non-stationary: Zt , the laboraugmenting technology process, and Υt , the investment-specic technology process. Zt grows at the gross rate πz t=Zt Zt−1 , while Υt grows at the gross rate πΥ t=Υt Υt−1 . The law of motion of all the disturbances is also described in the Appendix C. Let PX t be the nominal price of the variable X in period t . The consumption good is the numeraire and its price is Pt , which will be non-stationary. Let πt=Pt Pt−1 be the ination rate, Wt the nominal wage, RX t≡1 + rX t the gross interest rate in sector or variable X and rX t the net rate of interest. All other variables are expressed in real terms, unless otherwise stated. Finally, let e Xt denote the stationary version of X and e Xss its steady state. 5.2.1 Computing the stationary system manually Justiniano et al. (2011) state that ZtΥ α 1−α t is the balanced growth path for this economy. To show that this is indeed the case, we scale the optimality conditions of by this factor and show that they still hold after the transformation. From the nal and intermediate goods sectors' conditions we get: i) Final good production function: e At=e Qt. (28) ii) Pricing of nal good: e Pt=e PQ t. (29) iii) Intermediate goods production function: e Qt= (zL tLt)1−αe Ktα (30) iv) Wage and nominal rate denitions: f Wt≡Wt PtΥ α 1−α tZt ;e RK,t ≡ΥtRK,t Pt (31) 15 v) Labor demand function: Lt= (1 −α)g MCt f Wte Qt. (32) vi) Capital demand function: e Kt=αg MCt e RK te Qt. (33) vii) Using PQ t=PQ t(n) and πQ t≡PQ t PQ t−1 , the cost of adjusting prices is γP Q,t =ϕP Q 2"πQ t πQ t−1 −1#2 . (34) viii) The Euler equation for intermediate rms is: e Qt−θH te Qt+g MCtθH te Qt e PQ t −ϕP Q "πQ t πQ t−1 −1#πQ t πQ t−1e Qt +Et(ΛϕP Q "πQ t+1 πQ t −1#(πQ t+1)2 πQ te Qt+1 πΥ t+1α 1−απz t+1)= 0. (35) The rst order conditions for the investment and capital producer give: ix) Production function for investment goods: e It=e YI t. (36) x) Investment price: e PI t= 1. (37) xi) Production function for capital goods: e Knew t=zI,t 1−S e It e It−1!!e It. (38) xii) Optimality condition for investments: "Λ1 πΥ t+1 e PK t+1zI,t+1S′ 2 e It+1 e It!e It+1#− e PK tzI,t "1−S e It e It−1!−S′ 1 e It e It−1!e It#−e PI t= 0 (39) xiii) Investment adjustment costs: S e It e It−1!=ϕI1 2"e It e It−1πΥ t1 1−απz t−πΥ ss1 1−απz ss#2 . (40) 16 and thus S′ 1 e It e It−1!e It=ϕI1he It e It−1πΥ t1 1−απz t−πΥ ss1 1−απz ssi e It e It−1πΥ t1 1−απz t . (41) S′ 2 e It e It−1!e It=− ϕI1he It e It−1πΥ t1 1−απz t−πΥ ss1 1−απz ssi e It e It−1πΥ t1 1−απz t2. (42) The rst order conditions for the households and the entrepreneurs give: Xiv) Marginal utility of consumption: Υ α 1−α tZtu′(Ct)≡eu′e Ct=zu t   e Ct−e Ct−1bc (πΥ t)α 1−απz t 1−bc (πΥ ss)α 1−απz ss     −1 . (43) xv) Marginal disutility of labor: v′(Lt) = Lt(j)−blLt−1 1−blζ . (44) xvi) The stochastic discount factor: Λt,t+i≡βiρt+1eu′e Ct+i ρteu′e Ct1 πt+i 1 πΥ t+iα 1−απz t+i , (45) Then the one-period ahead stochastic discount factor is then Λt,t+1 ≡β ρt+1eu′e Ct+1 ρteu′e Ct1 πt+1 1 πΥ t+1α 1−απz t+1 . (46) xvii) Optimality condition for bond holdings: Et[Λt,t+1]Rt= 1, (47) xviii) Labor supply equation: f Wt=ψtMRS(Lt,e Ct)    [ψt−1] [1 −γt] +ϕWπW t πW t−1πW t πW t−1 −1 −EhΛϕWLt+1 Lt πW t+1 πW tπW t+1 πW t −1i    −1 . (48) where MRS(Lt, Ct) = Υ α 1−α tZtv′(Lt) eu′e Ct= Υ α 1−α tZtMRS(Lt,e Ct). (49) 17 xix) Wage adjustment equation: γt≡ϕW 2"πW t πW t−1 −1#2 . (50) xx) Optimality condition for capital: e PK t=Et Λπt+1   1 πΥ t+1 e PK t+1 (1 −δ) +1 πΥ t+1 e RK,t+1Ut+1 −1 πΥ t+1 eγ(Ut+1)  . (51) xxi) Optimality condition for capital utilization: e RK,t =eγ′(Ut), (52) where we have used that eγ(Ut)≡Υtγ(Ut) = e RK ss ϕuheϕu(Ut−1) −1i, (53) eγ′(Ut)≡Υtγ′(Ut) = e RK,sseϕu(Ut−1). (54) xxii) Eective capital denition: e Kt=Ute Kt−1 πΥ t1 1−απz t . (55) xxiii) Physical capital accumulation equation: e Kt=(1 −δ)e Kt−1 πΥ t1 1−απz t +e Knew t. (56) xxiv) Final goods market clearing: e At=e Ct+e YI t. (57) xxv) Final goods ination: πt≡Pt Pt−1 . (58) xxvi) Intermediate goods ination: πQ t=πte PQ t e PQ t−1 . (59) xxvii) Wage ination: πW t=f Wt f Wt−1 πtπz tπΥ tα 1−α. (60) xxviii) Output denition: e YNAT,t =e At. (61) 18 xxix) Output growth denition: △e YNAT,t =e YNAT,t e YNAT,t−1 πz tπΥ tα 1−α. (62) A le with these equations and a set of instructions used to solve the system with the NB toolbox are in the Appendix D. 5.2.2 Computing the stationary system automatically The steps are the same as in the RBC example. First, we write the equations of the non-stationary model to a le and save it with name jpt_non_stationary.nb. This le is reproduced in the Appendix E together with the commands needed to make the model stationary and compute impulse responses. The resulting le with the stationary model is also reproduced in the Appendix F. To check for the accuracy of the automated solution, we compare impulse responses to an investment-specic technology shock in the stationary version of the model computed manually and automatically and plot them in gure 2. The two responses are identical also in this case. An investment-specic shock permanently raises consumption, investment, and output although the level of consumption and level of output temporarily move in dierent directions. As an investment unit is permanently transformed into more capital, investment, consumption and output as a deviation from the steady state temporarily fall and converge slowly. The responses of stationary variables are persistent because price adjustment costs are important. Hours also fall but temporarily while price ination and wage ination initially fall and then rebound as demand for transitory investment builds up. Finally, given the Taylor rule, the dynamics of the interest rate track the dynamics of ination along the adjustment path. It is instructive to analyze the Gi functions produced. They are characterized by 18 equations: A△ t=Q△ t, (63) T△ t−α¯ K△ t= (1 −α)Z△ t, (64) ¯ K△ t+RK△ t=T△ t, (65) W P△ t =T△ t, (66) I△ t−Υ△ t=YI△ t, (67) PI P△ t =−Υ△ t, (68) Knew△ t=I△ t, (69) 19 Consumption gap 5 25 45 65 85 0.95 1 1.05 1.1 0.95 1 1.05 1.1 Inflation 5 25 45 65 85 1 1.01 1.02 1 1.01 1.02 Wage inflation 5 25 45 65 85 1.01 1.02 1.03 1.01 1.02 1.03 Investment gap 5 25 45 65 85 0.05 0.1 0.15 0.2 0.05 0.1 0.15 0.2 Capital gap 5 25 45 65 85 1 2 3 4 5 1 2 3 4 5 Hours worked 5 25 45 65 85 0.85 0.9 0.95 0.85 0.9 0.95 Output gap 5 25 45 65 85 1 1.1 1.2 1.3 1 1.1 1.2 1.3 Money market rate 5 25 45 65 85 1 1.01 1.02 1 1.01 1.02 Consumption 5 25 45 65 85 1 1.05 1.1 1.15 1 1.05 1.1 1.15 Investment 5 25 45 65 85 1 2 3 1 2 3Capital 5 25 45 65 85 1 1.5 2 2.5 1 1.5 2 2.5 Output 5 25 45 65 85 0.9 1 1.1 1.2 0.9 1 1.1 1.2 AUTOMATIC STATIONARY MANUAL STATIONARY Steady state Figure 2: Responses of variables to an temporary shock to investment-specic technology. The green line gives the steady state/inital value. PK P△ t =PI P△ t , (70) C△ t−u′(C)△ t= 0, (71) W P△ t =v′(L) u′(C)△ t , (72) v′(L) u′(C)△ t =−u′(C)△ t, (73) PK P△ t =RK△ t, (74) PK P△ t =γ(U)△ t, (75) 20 RK△ t=γ′(U)△ t, (76) K△ t=Knew△ t, (77) ¯ K△ t=K△ t, (78) A△ t=C△ t, (79) Y△ t=A△ t. (80) (63)-(64) come from the nal good and the intermediate production functions, (65)-(66) from the intermediate producer conditions, (67)-(68) from the investment good production function and optimization, (69)-(70) from capital production and optimization, (71)-(78) from the consumer optimality conditions, while (79)-(80) are obtained from market clearing and the resource constraint. Recall that these equations characterize the trending variables. Note that, while capital, its real rate and the real wage are on the list, hours are not. Also, since the relative price of capital is trending, the marginal rate of substitution between leisure and consumption is also trending. We contrast these facts next with those implied by a model with unit root in preferences. 5.2.3 Changing the location of the unit root To investigate the consequences of assuming unit roots in dierent disturbances, we let the preference shock, rather than the investment-specic technology shock, be a unit root process. In practice, this requires removing Υt from the list of processes generating non-stationary dynamics, and adding Xt to it, where Xt is the preference shock, entering the utility as follows u(Ct(j)) = Xtzu t1−bc πz ss ln  Ct(j)−bcCt−1 1−bc πz ss  , (81) The model le is in the Appendix G. The automatically stationarized le in the Appendix I. Clearly, when dierent disturbances feature unit root, the Gi function will generally change as dierent variables display trending behavior. Thus, a comparison of the Gi functions allows a researcher to understand the specication that is more likely to be compatible with the data. When the preference shock features a unit root equations (63), (69), (72), (77), (78), (79) and (80) are unchanged. Equations (67),(68), (70), (74)-(76) are no longer relevant. From the intermediate production function, we now have: (α−1)L△ t+T△ t−α¯ K△ t= (1 −α)Z△ t, (82) while equations (65), (66) and (71) become ¯ K△ t+RK△ t=T△ t, (83) 21 Table 1 Balanced growth, selected variables. Model with a unit root Model with a unit root in investment technology in preferences C 1.0030 1.0029 I 1.0090 1.0029 K 1.0090 1.0029 L 1.0000 1.0011 Y 1.0030 1.0029 PI P 0.9941 1.0000 W 1.0030 1.0018 L△ t+W P△ t =T△ t, (84) C△ t−u′(C)△ t=X△ t. (85) Compared with the version with investment specic unit root, the return on capital is now stationary while hours worked are trending to allow for the permanent increase in preferences for consumption. As a consequence, the marginal utility of labor is also trending and the marginal rate of substitution is also aected because with trending hours, real wages do not grow at same rate as output. ζL△ t=v′(L)△ t, (86) v′(L) u′(C)△ t =v′(L)△ t−u′(C)△ t. (87) Finally, the market clearing condition imply an additional restriction: A△ t=YI△ t. (88) We report the growth rates of few variables of the model in table 1. In the economy with investmentspecic technology shocks, investment and capital increase faster than output over time, as the real price of investment decreases over time, while this is not the case in the economy with unit roots in preferences. Thus, this set of facts can be used to select which specication is more in line with the data. Alternatively, suppose that the labor share is of interest. If, in the data, it is trending and hours are stationary, which happens in US, France, Germany and Canada (see e.g. Canova and Matthes (2021)) then a model with unit root in investment-specic shocks should be preferred, as it implies trending output and stationary hours. If, on the other hand, the labor share and hours are both trending, as it happens, for example, in the UK, then a model with unit root in preferences is more appealing, as both output and hours display a 22 References Blanchard, O. J. and Kahn, C. M. (1980) The Solution of Linear Dierence Models under Rational Expectations, Econometrica , 48 , 13051311. Canova, F. (2014) Bridging dsge models and the raw data, Journal of Monetary Economics , 67 , 115. Canova, F. and Matthes, C. (2021) A composite likelihood approach fro dynamic structural models, Economic Journal , 131 , 24472477. Comin, D. and Gertler, M. (2006) Medium term business cycles, American Economic Review , 96 , 523551. Dynare (2021) https://www.dynare.org/ , Dynare, CEPREMAP. GPMN (2020) https://iris.igpmn.org/ , Iris macroeconomic modeling toolbox, Global Projection Model Network. Justiniano, A., Primiceri, G. and Tambalotti, A. (2011) Investment shocks and the relative price of investment, Review of Economic Dynamics , 14 , 101121. Klein, P. (2000) Using the generalized shur form to solve a multivariate linear rational expectation model, Journal of Economic Dynamics and Control , 24 , 14051423. Lafourcade, P. and de Wind, J. (2012) Taking Trends Seriously in DSGE Models: An Application to the Dutch Economy, DNB Working Papers 345, Netherlands Central Bank, Research Department. Maih, J. (2021) https://github.com/jmaih/RISE_toolbox , Rise, Norges Bank. Paulsen, K. S. (2021) https://github.com/Coksp1/NBTOOLBOX , toolbox, Norges Bank. Sims, C. (2010) Gensys, manuscript, Princeton University. 29 A Description of how the algorithm is implemented The algorithm utilizes dierent classes existing in the NB toolbox: 1. The model le is parsed, and written into a MATLAB function handle. 2. The process of nding the restrictions imposed on the balanced growth path is performed utilizing the nb_bgrowth class. First, endogenous and exogenous variables, unit root process and parameters are transformed to a set of nb_bgrowth objects which are passed to the MATLAB function handle representing the equations of the model. The result is the set of restrictions holding along the balanced growth path. 3. The initial solution for the balanced growth path is obtained by transforming the equations found in step 2 to a linear system. From this solution it is possible to identify the R non-stationary variables. 4. The identication of the R the Gi functions is done with the function nb_findBasis . To get the symbolic representation of the Gi 's, we rst construct a function handle representing these equations, and then utilize the nb_term class and subclasses to transform it into symbolic equations. 5. The stationary representation of the model is constructed using nb_stTerm . This requires transforming the endogenous and exogenous variables and unit root variables to class nb_stTerm objects using the solution for the balanced growth. The parameters are transformed to nb_stParam objects. Then these objects are passed to the MATLAB function handle representing the non-stationary equations. The result is the stationary representation written in the model le syntax. 6. With the stationary representation, the model can be solved in a standard way. The examples presented in the paper utilizes the myAD 2 , or nb_mySD classes. Again, endogenous and exogenous variables, unit root process and parameters are transformed to a set of objects which are passed to the MATLAB function handle representing the stationary equations of the model. The result is a set of matrices representing the linearized system, which is solved by Klein (2000) algorithm. Step 1 uses the nb_dsge class constructor; steps 2-4 the nb_dsge.solveBalancedGrowthPath method 3 , while step 5 employs the nb_dsge.stationarize method. The nal step is performed by the nb_dsge.solve method. The following command allow one to perform model-based analysis once the solution is found: IRFs Produces impulse response functions using the method nb_dsge.irf . Theoretical moments Produce the theoretical moments of the model using the method nb_dsge.theoreticalMoments . 2 Original code is written by SeHyoun Ahn, and can be found here https://github.com/sehyoun/MATLABAutoDiff 3 A method is a function that is related to a class. For example the nb_dsge.solveBalancedGrowthPath method is a function that act on an object of class nb_dsge . 30 Filter If nb_dsge is assigned to data, using the nb_dsge.set method, you can run the Kalman lter using the nb_dsge.filter method. See also the nb_dsge.getFiltered method. Forecast If the model is ltered, forecast can be produced using the nb_dsge.forecast method. See also nb_dsge.getForecast and nb_dsge.plotForecast methods. B The details of the RBC model The production function is Yt= (AtLt)γK1−γ t−1, (116) where Yt is output, Lt= 1 is labor, Kt is capital, At is the labor augmenting productivity, and γ is the labor share. Let It be investment and δ is the depreciation rate. Capital accumulates according to Kt= (1 −δ)Kt−1+It, (117) A representative household maximizes the discounted utility E"∞ X s=0 βsU(Ct+s)#, (118) where β is the discount factor, subject to (116), (117) and the constraint Bt+Yt=Ct+It(1 + rt−1)Bt−1+Tt, (119) where Bt are one period bond issued by the government. The government collects lump sum taxes, Tt , from the households. The optimality condition is E[Ct+1] Ct =β(1 + rt), (120) where rt , the real interest rate, is given by rt= (1 −γ)E[Yt+1] Kt −δ. (121) Markets must also clear Yt=Ct+It. (122) (120)-(121) together with (116), (117) , (122), determine the equilibrium. The stationary model le 31 endogenous y r l k i dA c D_Z_y D_Z_k D_Z_i D_Z_c D_Z_A exogenous u parameters std_u lambda gamma g delta beta model (c(+1)*D_Z_c(+1))/c = beta*(1+r); l-1; y = ((1*l)^gamma)*((k(-1)*D_Z_k^-1)^(1-gamma)); ((1-gamma)*(y(+1)*D_Z_y(+1)))/k = r+delta; k = (1-delta)*(k(-1)*D_Z_k^-1)+i; y = c+i; dA = 1/(1*D_Z_A^-1); 1/(1*D_Z_A^-1) = ((g^(1-lambda))*(((1*D_Z_A^-1)/(1*D_Z_A^-1*D_Z_A(-1)^-1))^lambda))*(exp(std_u*u)); gamma*log(D_Z_k)+log(D_Z_y)-gamma*log(D_Z_A)-log(D_Z_k); log(D_Z_y)-log(D_Z_k); log(D_Z_k)-log(D_Z_i); log(D_Z_y)-log(D_Z_c); reporting A=exp(cumsum(log(dA) - log(steady_state(dA)))); C=c*A; I=i*A; K=k*A; Y=y*A; c_log_dev=exp(log(c) - log(steady_state(c))); i_log_dev=exp(log(i) - log(steady_state(i))); k_log_dev=exp(log(k) - log(steady_state(k))); r_dev=r - steady_state(r); y_log_dev=exp(log(y) - log(steady_state(y))); C The details of the New Keynesian model Final good sector The nal goods sector combines intermediate goods, Qt(n) into a nal good At sold at a price Pt . The production function is: Qt=Z1 0 Qt(n)1−1 θH tdnθH t θH t−1. (123) where θH t is the elasticity of substitution between dierent intermediate goods, and it assumed to follow an AR process. The optimal combination of Qt(n) is found by cost minimization 32 min {Qt(n)}Z1 0 PQ t(n)Qt(n)dn, subject to (123). The resulting optimality condition is: Qt(n) = PQ t(n) PQ t!−θH t Qt, (124) where PQ t=Z1 0 PQ t(n)1−θH tdn1 1−θH t. (125) Final good rms maximize the prot function PAt−PQ tQt, (126) subject to At=Qt. (127) The optimality condition is Pt=PQ t. (128) Intermediate goods sector Intermediate goods rms use eective capital and labor to produce a good which is sold under monopolistic competition to the nal goods sector. The intermediate rm n has the production function Qt(n)=(ZtzL tLt(n))1−αKt(n)α, (129) where α∈[0,1] is the capital share, Lt(n) and Kt(n) denote hours and eective capital of rm n in period t . There are two shocks to productivity: Zt , a permanent labor-augmenting technology process, growing at the gross rate πz t , and zL t , a temporary shock to productivity (or labor utilization). Total labor input to rm n aggregate labor inputs from all households j : Lt(n) =   1 Z0 Lt(n, j)1−1 ψtdj  ψt ψt−1 , (130) 33 where ψt is the elasticity of substitution between dierentiated labor. Let Wt be the wage rate and RK t the rental rate. Minimizing total factor outlays: WtLt(n) + RK tKt(n), subject to 129 leads to the following rst order conditions Lt= (1 −α)MCt Wt Qt(n), (131) Kt=αMCt RK t Qt(n), (132) where MCt is marginal cost. In symmetric equilibrium all n rms make the same decision, so Lt=Lt(n) , Kt=Kt(n) and Qt=Qt(n) . Firm n will minimize labor costs subject to the optimal level of labor input: min lI,t(n,j)Z1 0 Lt(n, j)Wt(j)dj, (133) subject to   1 Z0 Lt(n, j)1−1 ψtdj  ψt ψt−1 =Lt(n). (134) Thus the conditional labor demand functions facing household j is Lt(j) = Wt(j) Wt−ψt Lt. (135) Firms in the intermediate sector sell their goods monopolistically. Each rm n charges the PQ t(n) . Prots paid out as dividends to households) are: Πt(n) = PQ t(n)Qt(n)−WtLt(n)−RK t(n)Kt(n). (136) The costs of adjusting prices are γP Q,t(n)≡ϕP Q 2"PQ t(n)/PQ t−1(n) PQ t−1/PQ t−2 −1#2 , (137) where PQ t=hR1 0PQ t(n)1−θH tdni1 1−θH t . The costs of changing prices is governed by the parameter ϕP Q . These costs are assumed to be intangible. The intermediate rm n faces the demand function: Qt(n) = PQ t(n) PQ t−θH t Qt Optimal price setting for rm n requires maximizing 4 : 4 To make the expression easier to work with, the costs of adjusting prices are linear in PQ tQt and not PQ t(n)Qt(n) 34 Πs= ∞ X t=s ∆s,t      PQ t(n)PQ t(n) PQ t−θH t Qt−MCt"PQ t(n) PQ t−θH t Qt# −γP Q,t(n)PQ tQt     . In symmetric equilibrium, all rms will behave the same, and thus: Qt−θH tQt+MCtθH t Qt PQ t −ϕP Q "PQ t/PQ t−1 PQ t−1/PQ t−2 −1#PQ t/PQ t−1 PQ t−1/PQ t−2 Qt +Et(∆ϕP Q "PQ t+1/PQ t PQ t/PQ t−1 −1#(PQ t+1/PQ t)2 PQ t/PQ t−1 Qt+1)= 0. (138) Investment producer Perfectly competitive rms purchase YI t units from the nal good producers, and produce investment goods It according to the production function It= ΥtYI t. (139) Υt represents investment-specic technological (IST) progress, which is specied later, and It is the investment good sold to capital producers. The objective of investment producers is to maximize the prot function PI tIt−PtYI t, (140) subject to (139). The optimal condition is: ΥtPI t=Pt. (141) Capital producer Capital goods, Knew t , are produced by a sector which purchases investment goods It and transform them into capital, sold to households at price PK t . The production function is: Knew t=zI,t 1−SIt It−1It. (142) S captures the presence of adjustment costs in investment. zI,t is the marginal eciency of investment shock process. The objective of capital producers is to maximize the expected discounted future prots: Es ∞ X t=s ∆s,t PK tKnew t−PI tIt, (143) 35 subject to (139). Optimality implies: E∆PK t+1zI,t+1S′ 2It+1 ItIt+1+PK tzI,t 1−SIt It−1−S′ 1It It−1It−PI t= 0, (144) where we assume that SIt It−1=ϕI1 2It It−1 −πΥ ss1 1−απz ss2 , (145) Thus S′ 1It It−1=ϕI1It It−1 −πΥ ss1 1−απz ss1 It−1 , (146) S′ 2It It−1=−ϕI1It It−1 −πΥ ss1 1−απz ssIt I2 t−1 , (147) Households Each household supplies a dierentiated labor input to intermediate rms and set wages under the assumption of monopolistic competition. Households obtain utility from consumption and leisure. Preferences are additively separable. Lifetime expected utility of household j at time s is Us(j) = Es ∞ X t=s βt−sρt[u(Ct(j)) −v(Lt(j))] , (148) where β is the discount factor, ρt is a discount factor shock, Ct denotes consumption, and Lt labor. The in-period utility function is u(Ct(j)) = zu t 1−bc (πΥ ss)α 1−απz ss !ln     Ct(j)−bcCt−1 1−bc (πΥ ss)α 1−απz ss    , (149) ν(Lt(j)) = 1−bl 1 + ζLt(j)−blLt−1 1−bl1+ζ , (150) bc governs habit persistence, πz ss denotes the steady-state labor augmenting technology growth rate and πΥ ss the steady-state investment-specic technology growth rate. The degree of disutility of supplying labor is captured by ζ > 0 , the inverse Frisch elasticity. As (149) indicates, we assume a log in-period utility function for consumption. This insure the existence of a balanced growth path. The household's budget constraint is 36 PtCt(j) + PK tKnew t(j) + PtBt(j) + Ptγ(Ut)Kt−1 +Wt(j)Lt(j)γt(j) + Tt(j) ≤Rt−1Pt−1Bt−1(j) + Wt(j)Lt(j) + RK tUtKt−1 (151) +Qt(j)+Πt(j) + DIVt(j) , where Pt is the price level of nal goods, Rt is the gross interest rate, Bt(j) is real household's borrowing, Wt(j) is the nominal wage rate set by household j , γt(j) is intangible wage adjustment costs (dened in (156)), Lt(j) is the total hours worked, and DIVt(j) , Πt(j) and Tt(j) are is the net cash ow from household's j portfolio of state contingent securities, prots (in nominal terms) paid out and taxes paid, respectively. Households rent out utilized capital UtKt−1 at the rate RK t , where Ut is the utilization rate and Kt−1 the capital stock at the end of last period. We assume intangibility of the utilization cost. Knew t(j) is new capital bought from the capital producers at price PK t and γ(Ut)Kt−1 is a cost specied and: γ(Ut) = RK ss Pssϕuheϕu(Ut−1) −1i, (152) where ϕu governs the cost of adjusting utilization. Time t utilized capital iis Kt=UtKt−1. (153) The physical capital accumulation equation is Kt= (1 −δ)Kt−1+Knew t. (154) Household j faces the following labor demand curve: Lt(j) = Wt(j) Wt−ψt Lt , (155) where Wt is the wage rate. We further assume that there is sluggish wage adjustment due to resource costs measured in terms of the total wage bill. We deviate from JPT, and assume that wage adjustments costs are γt(j)≡ϕW 2Wt(j)/Wt−1(j) Wt−1/Wt−2 −12 . (156) Thus, costs are related to changes in wage ination relative to the past observed rate for households. ϕW>0 determines how costly it is to change the wage ination rate. Maximizing (148) subject to the budget constraint leads to the following optimal conditions in symmetric equilibrium 37 Etβρt+1u′(Ct+1) ρtu′(Ct) Pt Pt+1 Rt= 1 (157) v′(Lt) u′(Ct)ψt Pt Wt =(ψt−1) (1 −γt) + ϕWWt/Wt−1 Wt−1/Wt−2 −1Wt/Wt−1 Wt−1/Wt−2 −Etβρt+1u′(Ct+1) ρtu′(Ct) Pt Pt+1 Lt+1 Lt ϕWWt+1/Wt Wt/Wt−1 −1(Wt+1/Wt)2 Wt/Wt−1, (158) ρtu′(Ct)PK t Pt =Et"βρt+1u′(Ct+1) PK t+1 Pt+1 (1 −δ) + RK,t+1 Pt+1 Ut+1 −γ(Ut+1)!#, (159) RK,t Pt =γ′(Ut)RK,ss Pss eϕu(Ut−1). (160) Final good market clearing and denitions e At=e Ct+e YI t. (161) Ination πt≡Pt Pt−1 . (162) Intermediate goods ination πQ t=PQ t PQ t−1 . (163) Investment price ination πI t=PI t PI t−1 . (164) Capital price ination πK t=PK t PK t−1 . (165) Wage ination πW t=Wt Wt−1 . (166) Output YNAT,t =At. (167) 38 DPQ_Y_NW = (NAT_Y_NW/NAT_Y_NW(-1))*DZT_NW*DUT_NW^(ALPHA_NW/(1 - ALPHA_NW)); % Intermidate good inflation (DPQ_PQ_NW) DPQ_PQ_NW = DPQ_P_NW*REAL_PQ_NW/REAL_PQ_NW(-1); % Real investment inflation (DPQ_REAL_PI_NW) DPQ_REAL_PI_NW = (REAL_PI_NW/REAL_PI_NW(-1))/DUT_NW; % Real wage inflation (DPQ_REAL_W_NW) DPQ_REAL_W_NW = (REAL_W_NW/REAL_W_NW(-1))*DZT_NW*DUT_NW^(ALPHA_NW/(1 - ALPHA_NW)); % Wage inflation (DPQ_W_NW) DPQ_W_NW = REAL_W_NW/REAL_W_NW(-1)*DPQ_P_NW*DZT_NW*DUT_NW^(ALPHA_NW/(1 - ALPHA_NW)); % 7) Monetary policy rule (DPQ_P_NW) RN3M_NW = RN3M_NW(-1)^OMEGA_R_NW*(steady_state(RN3M_NW) *(DPQ_P_NW/steady_state(DPQ_P_NW))^OMEGA_P_NW *(NAT_Y_NW/steady_state(NAT_Y_NW))^OMEGA_Y_NW *(DPQ_Y_NW/steady_state(DPQ_Y_NW))^OMEGA_DPQ_Y_NW) ^(1-OMEGA_R_NW)*(Z_RN3M_NW/steady_state(Z_RN3M_NW)); [static] DPQ_P_NW = DPQ_P_NW_SS; % 8) Shock processes % Permanent labor-augmenting technology process (DZT_NW) log(DZT_NW) = (1-LAMBDA_DZT_NW)*log(DZT_NW_SS) + LAMBDA_DZT_NW*log(DZT_NW(-1)) + E_DZT_NW*std_E_DZT_NW; % Investment-specific technological progress (DUT_NW) log(DUT_NW) = (1-LAMBDA_DUT_NW)*log(DUT_NW_SS) + LAMBDA_DUT_NW*log(DUT_NW(-1)) + E_DUT_NW*std_E_DUT_NW; % Competition in the labor market shock process (PSI_NW) log(PSI_NW) = (1-LAMBDA_PSI_NW)*log(PSI_NW_SS) + LAMBDA_PSI_NW*log(PSI_NW(-1)) + E_PSI_NW*std_E_PSI_NW ; % Discount factor shock process (RHO_NW) log(RHO_NW) = (1-LAMBDA_RHO_NW)*log(RHO_NW_SS) + LAMBDA_RHO_NW*log(RHO_NW(-1)) + E_RHO_NW*std_E_RHO_NW ; % Price markup shock process (THETAH_NW) log(THETAH_NW) = (1-LAMBDA_THETAH_NW)*log(THETAH_NW_SS) + LAMBDA_THETAH_NW*log(THETAH_NW(-1)) + E_THETAH_NW*std_E_THETAH_NW ; % Marginal efficiency of investment shock process (Z_I_NW) log(Z_I_NW) = (1-LAMBDA_I_NW)*log(Z_I_NW_SS) + LAMBDA_I_NW*log(Z_I_NW(-1)) + E_I_NW*std_E_I_NW ; % Temporary labor augmenting technology shock process (Z_L_NW) log(Z_L_NW) = (1-LAMBDA_L_NW)*log(Z_L_NW_SS) + LAMBDA_L_NW*log(Z_L_NW(-1)) + E_L_NW*std_E_L_NW ; 45 % Consumption preference shock process (Z_U_NW) log(Z_U_NW) = (1-LAMBDA_U_NW)*log(Z_U_NW_SS) + LAMBDA_U_NW*log(Z_U_NW(-1)) + E_U_NW*std_E_U_NW ; % Monetary policy shock process (Z_RN3M_NW) log(Z_RN3M_NW) = (1-LAMBDA_RN3M_NW)*log(Z_RN3M_NW_SS) + LAMBDA_RN3M_NW*log(Z_RN3M_NW(-1)) + E_RN3M_NW*std_E_RN3M_NW; % Construct reported variables that can be asked for in IRFs, but % are not part of the model. Here you can use all MATLAB functions % that act on a double vector, and returns an output with the same % size as the input. reporting C_NW_LEVEL = exp(cumsum(log(DPQ_C_NW) - log(steady_state(DPQ_C_NW)))); I_NW_LEVEL = exp(cumsum(log(DPQ_I_NW) - log(steady_state(DPQ_I_NW)))); K_NW_LEVEL = exp(cumsum(log(DPQ_K_NW) - log(steady_state(DPQ_K_NW)))); Y_NW_LEVEL = exp(cumsum(log(DPQ_Y_NW) - log(steady_state(DPQ_Y_NW)))); To parse and solve this model you can use the following instructions: % Read the stationary model modelS = nb_dsge('nb_file','jpt.nb'); % Give it a name modelS = set(modelS,'name','Stationarize JPT manually'); % Set the parameters param = struct(); param.ALPHA_NW = 0.167; param.BC_NW = 0.859; param.BETA_NW = 100/(0.134+100); param.BL_NW = 0; param.DELTA_NW = 0.025; param.DPQ_P_NW_SS = (0.702 + 100)/100; param.DUT_NW_SS = 1 + (0.597/100); param.DZT_NW_SS = 1 + (0.303 − (param.ALPHA_NW/... (1 − param.ALPHA_NW))*0.597)/100; param.LAMBDA_DUT_NW = 0.156; param.LAMBDA_DZT_NW = 0.286; param.LAMBDA_I_NW = 0.772; param.LAMBDA_L_NW = 0; param.LAMBDA_PSI_NW = 0.967; param.LAMBDA_RHO_NW = 0.590; param.LAMBDA_RN3M_NW = 0; param.LAMBDA_U_NW = 0; param.LAMBDA_THETAH_NW = 0.971; param.OMEGA_DPQ_Y_NW = 0.208; param.OMEGA_P_NW = 1.709; param.OMEGA_R_NW = 0.858; param.OMEGA_Y_NW = 0.051; param.PHI_PQ_NW = 0.2; param.PHI_I1_NW = 2.657; 46 param.PHI_W_NW = 1.0080; param.PHI_U_NW = 5.434; param.PSI_NW_SS = 1.135/(1.135 − 1); param.RHO_NW_SS = 1; param.THETAH_NW_SS = 1.171/(1.171 − 1); param.ZETA_NW = 4.444; param.Z_I_NW_SS = 1; param.Z_L_NW_SS = 1; param.Z_U_NW_SS = 1; param.Z_RN3M_NW_SS = 1; param.std_E_DUT_NW = 0.630; param.std_E_DZT_NW = 0.933; param.std_E_I_NW = 5.103; param.std_E_L_NW = 0; param.std_E_PSI_NW = 0.310; param.std_E_RHO_NW = 0.036; param.std_E_RN3M_NW = 0.210; param.std_E_THETAH_NW = 0.219; param.std_E_U_NW = 0; modelNS = assignParameters(modelNS,param); % Solve steady state numerically ssInit = struct(... 'GAMMA_W_NW',0,... 'PSI_NW',param.PSI_NW_SS,... 'S_NW',0,... 'S_PRIME1_NW',0,... 'S_PRIME2_NW',0,... 'GAMMA_U_NW',0,... 'THETAH_NW',param.THETAH_NW_SS); modelS = checkSteadyState(modelS,... 'solver','fsolve',... 'steady_state_default', @ones,... 'steady_state_init', ssInit,... 'steady_state_solve', true); ss = getSteadyState(modelS) % Solve stationary model modelS = solve(modelS); 47 E The non-stationary model le endogenous A_NW % Final goods production C_NW % Consumption DPQ_C_NW % Consumption growth DPQ_I_NW % Investment growth DPQ_K_NW % Capital growth DPQ_P_NW % Inflation DPQ_PQ_NW % Intermidate goods inflation DPQ_REAL_PI_NW % Real investment inflation DPQ_REAL_W_NW % Real wage inflation DPQ_W_NW % Wage inflation DPQ_Y_NW % Output growth DSA_NW % Stochastic discount factor GAMMAPRIME_U_NW % Marginal cost of utilizing the capital GAMMA_U_NW % Cost of utilize the capital GAMMA_W_NW % Wage adjustment cost I_NW % Investment K_NW % Capital KBAR_NW % Utilized capital KNEW_NW % Capital goods produced each period L_NW % Hours worked MC_NW % Marginal cost MRS_NW % Marginal rate of substitution NAT_Y_NW % Output PSI_NW % Competition in the labor market shock process Q_NW % Demand for intermediate goods REAL_PI_NW % Real price of investment REAL_PK_NW % Real price of capital REAL_PQ_NW % Real intermediate goods price REAL_W_NW % Real wage rate RHO_NW % Discount factor shock process RK_NW % Rental rate of capital RN3M_NW % Money market interest rate S_NW % Investment adjustment cost function S_PRIME1_NW % Derivative of the investment % adjustment cost function with respect to first input times I_NW S_PRIME2_NW % Derivative of the investment % adjustment cost function with respect to second input times I_NW T_NW % Intermediate goods production THETAH_NW % Price markup shock process U_NW % Utilization rate UPRIME_NW % Derivative of the utility function % of households wrt consumption VPRIME_NW % Derivative of the utility function % of households wrt labor Y_I_NW % Input to investment production Z_I_NW % Marginal efficiency of investment shock process Z_L_NW % Temporary labor augmenting technology shock process Z_RN3M_NW % Monetary policy shock process 48 Z_U_NW % Consumption preference shock process exogenous E_DUT_NW % Investment-specific technological innovation E_DZT_NW % Permanent labor-augmenting technology innovation E_I_NW % Marginal efficiency of investment innovation E_L_NW % Temporary labor augmenting technology innovation E_PSI_NW % Competition in the labor market innovation E_RHO_NW % Discount factor innovation E_RN3M_NW % Monetary policy innovation E_THETAH_NW % Price markup innovation E_U_NW % Consumption preference innovation parameters ALPHA_NW % Capital share BC_NW % Habit in consumption BETA_NW % Discount factor BL_NW % Habit in hours worked DELTA_NW % Depreciation rate DPQ_P_NW_SS % Steady-state inflation DUT_NW_SS % Steady-state growth rate in investment- % specific technology DZT_NW_SS % Steady-state growth rate in labor- % augmenting technology LAMBDA_DUT_NW % Shock persistent parameter % for the investment-specific technology shock LAMBDA_DZT_NW % Shock persistent parameter % for the labor-augmenting technology shock LAMBDA_I_NW % Shock persistent parameter % for the marginal efficiency of investment shock LAMBDA_L_NW % Shock persistent parameter % for the temporary labor augmenting technology shock LAMBDA_PSI_NW % Shock persistent parameter % for the competition in the labor market shock LAMBDA_RHO_NW % Shock persistent parameter % for the discount factor shock LAMBDA_RN3M_NW % Shock persistent parameter % for the monetary policy shock LAMBDA_U_NW % Shock persistent parameter % for the price markup shock LAMBDA_THETAH_NW % Shock persistent parameter % for the consumption preference shock OMEGA_Y_NW % Taylor rule coefficient on ouput gap. OMEGA_DPQ_Y_NW % Taylor rule coefficient on ouput growth gap. OMEGA_P_NW % Taylor rule coefficient on inflation gap. OMEGA_R_NW % Interest rate smoothing in the taylor rule PHI_PQ_NW % Intermediate goods price adjustment cost parameter PHI_I1_NW % Investment adjustment cost parameter PHI_W_NW % Wage adjustment cost parameter PHI_U_NW % Capital utilization cost parameter PSI_NW_SS % Steady-state elasticity of substitution % between differentiated labor 49 RHO_NW_SS % Discount factor shock process in steady state THETAH_NW_SS % Steady-state elasticity of % substitution between intermidate goods ZETA_NW % Inverse Frisch elasticity Z_I_NW_SS % Marginal efficiency of investment % shock in steady-state Z_L_NW_SS % Temporary labor augmenting technology % shock in steady-state Z_RN3M_NW_SS % Monetary policy shock in % steady-state Z_U_NW_SS % Consumption preference shock in % steady-state std_E_DUT_NW % Standard deviation of the % innovation to the Z_DUT_NW shock process std_E_DZT_NW % Standard deviation of the % innovation to the Z_DZT_NW shock process std_E_I_NW % Standard deviation of the % innovation to the Z_I_NW shock process std_E_L_NW % Standard deviation of the % innovation to the Z_L_NW shock process std_E_PSI_NW % Standard deviation of the % innovation to the PSI_NW shock process std_E_RHO_NW % Standard deviation of the % innovation to the RHO_NW shock process std_E_RN3M_NW % Standard deviation of the % innovation to the Z_RN3M_NW shock process std_E_THETAH_NW % Standard deviation of the % innovation to the THETAH_NW shock process std_E_U_NW % Standard deviation of the % innovation to the Z_U_NW shock process unitrootvars Z UPSILON model % 1) Final goods sector % Production function (A_NW) A_NW = Q_NW; % FOC (Q_NW) REAL_PQ_NW = 1; % 2) Intermediate goods sector % Intermidate production function (KBAR_NW) T_NW = (Z*Z_L_NW*L_NW)^(1-ALPHA_NW)*KBAR_NW^ALPHA_NW; % Optimality condition wrt utilized capital (I.e. demand function) (MC_NW) KBAR_NW = ALPHA_NW*(MC_NW/RK_NW)*T_NW; % Optimality condition wrt aggregated labor (I.e. demand function) (L_NW) 50 L_NW = (1 - ALPHA_NW)*(MC_NW/REAL_W_NW)*T_NW; % Pricing (REAL_PQ_NW) Q_NW - THETAH_NW*Q_NW + MC_NW*THETAH_NW*Q_NW/REAL_PQ_NW - 100*PHI_PQ_NW*(DPQ_PQ_NW/DPQ_PQ_NW(-1) - 1)*DPQ_PQ_NW/DPQ_PQ_NW(-1)*Q_NW + DSA_NW*100*PHI_PQ_NW*(DPQ_PQ_NW(+1)/DPQ_PQ_NW - 1) *DPQ_PQ_NW(+1)^2/DPQ_PQ_NW*Q_NW(+1) = 0; % 3) Investment producer % Investment production function (Y_I_NW) I_NW = UPSILON*Y_I_NW; % First order condition for investment (REAL_PI_NW) UPSILON*REAL_PI_NW = 1; % 4) Capital producer % Capital production function (KNEW_NW) KNEW_NW = Z_I_NW*(1 - S_NW)*I_NW; % Investment adjustment cost function (S_NW) S_NW = (PHI_I1_NW/2)*(I_NW/I_NW(-1) - bgp(I_NW))^2; % Derivative of the investment adjustment cost function with respect to first input % (Multiplied by I_NW!) (S_PRIME1_NW) S_PRIME1_NW = PHI_I1_NW*(I_NW/I_NW(-1) - bgp(I_NW))*I_NW/I_NW(-1); % Derivative of the investment adjustment cost function with respect to second input % (Multiplied by I_NW!) (S_PRIME2_NW) S_PRIME2_NW = -PHI_I1_NW*(I_NW/I_NW(-1) - bgp(I_NW))*(I_NW/I_NW(-1))^2; % Optimal capital investment (I_NW) 0 = DSA_NW*REAL_PK_NW(+1)*Z_I_NW(+1)*S_PRIME2_NW(+1) + REAL_PK_NW*Z_I_NW*(1 - S_NW - S_PRIME1_NW) - REAL_PI_NW; % 5) Households % Marignal utility of consumption (UPRIME_NW) UPRIME_NW = Z_U_NW*((C_NW - C_NW(-1)*BC_NW)/(1 - BC_NW/bgp(C_NW)))^(-1) ; % Marginal utility of labor labor (VPRIME_NW) VPRIME_NW = ((L_NW - BL_NW*L_NW(-1))/(1-BL_NW))^ZETA_NW ; % Stochastic discount factor (DSA_NW) DSA_NW = BETA_NW*RHO_NW(+1)*UPRIME_NW(+1)/(DPQ_P_NW(+1)*RHO_NW*UPRIME_NW); % Consumption euler equation (RN3M_NW) DSA_NW*RN3M_NW = 1; % Optimal wage setting (REAL_W_NW) REAL_W_NW = PSI_NW*MRS_NW/((PSI_NW-1)*(1-GAMMA_W_NW) + 1000*PHI_W_NW*DPQ_W_NW/DPQ_W_NW(-1)*(DPQ_W_NW/DPQ_W_NW(-1) - 1) - DSA_NW*DPQ_W_NW(+1)*L_NW(+1)/L_NW*1000*PHI_W_NW*DPQ_W_NW(+1) /DPQ_W_NW*(DPQ_W_NW(+1)/DPQ_W_NW - 1)); 51 % Wage adjusment cost (GAMMA_W_NW) GAMMA_W_NW = 1000*PHI_W_NW/2*(DPQ_W_NW/DPQ_W_NW(-1) - 1)^2 ; % Marginal rate of substitution between consumption and leisure (MRS_NW) MRS_NW = VPRIME_NW/UPRIME_NW; % 6) Entrepreneurs % optimality with respect to capital (REAL_PK_NW) REAL_PK_NW = DSA_NW*DPQ_P_NW(+1)*( REAL_PK_NW(+1)*(1 - DELTA_NW) + RK_NW(+1)*U_NW(+1) - GAMMA_U_NW(+1)); % optimality with respect to utilization (RK_NW) RK_NW = GAMMAPRIME_U_NW; % Cost of utilizing the capital (GAMMA_U_NW) GAMMA_U_NW = steady_state(RK_NW)/PHI_U_NW*(exp(PHI_U_NW*(U_NW - 1)) - 1) ; % Marginal cost of utilizing capital (GAMMAPRIME_U_NW) GAMMAPRIME_U_NW = steady_state(RK_NW)*exp(PHI_U_NW*(U_NW - 1)) ; % Capital accumulation (K_NW) K_NW = (1-DELTA_NW)*K_NW(-1) + KNEW_NW; % Capital utilization (U_NW) KBAR_NW = U_NW*K_NW(-1); % 7) Market clearing % Final good market clearing (C_NW) A_NW = C_NW + Y_I_NW; % Intermediate good market clearing (T_NW) T_NW = Q_NW; % 8) Definitions % Definition of natural output NAT_Y_NW = A_NW; % Consumption growth (DPQ_C_NW) DPQ_C_NW = C_NW/C_NW(-1); % Investment growth (DPQ_I_NW) DPQ_I_NW = I_NW/I_NW(-1); % Capital growth (DPQ_K_NW) DPQ_K_NW = K_NW/K_NW(-1); % Definition of natural output growth (DPQ_A_NW) DPQ_Y_NW = NAT_Y_NW/NAT_Y_NW(-1); % Definition of intermidate good inflation (DPQ_PQ_NW) DPQ_PQ_NW = DPQ_P_NW*REAL_PQ_NW/REAL_PQ_NW(-1); 52 % Real investment inflation (DPQ_REAL_PI_NW) DPQ_REAL_PI_NW = REAL_PI_NW/REAL_PI_NW(-1); % Real wage inflation (DPQ_REAL_W_NW) DPQ_REAL_W_NW = REAL_W_NW/REAL_W_NW(-1); % Definition of wage inflation (DPQ_W_NW) DPQ_W_NW = REAL_W_NW/REAL_W_NW(-1)*DPQ_P_NW; % 9) Taylor rule (DPQ_P_NW) RN3M_NW = RN3M_NW(-1)^OMEGA_R_NW*(steady_state(RN3M_NW) *(DPQ_P_NW/steady_state(DPQ_P_NW))^OMEGA_P_NW *(NAT_Y_NW/steady_state(NAT_Y_NW))^OMEGA_Y_NW *(DPQ_Y_NW/steady_state(DPQ_Y_NW))^OMEGA_DPQ_Y_NW) ^(1-OMEGA_R_NW)*(Z_RN3M_NW/steady_state(Z_RN3M_NW)); [static] DPQ_P_NW = DPQ_P_NW_SS; % 10) Shock processes % Labor market competition shock (PSI_NW) log(PSI_NW) = (1-LAMBDA_PSI_NW)*log(PSI_NW_SS) + LAMBDA_PSI_NW*log(PSI_NW(-1)) + E_PSI_NW*std_E_PSI_NW ; % Discount factor shock (RHO_NW) log(RHO_NW) = (1-LAMBDA_RHO_NW)*log(RHO_NW_SS) + LAMBDA_RHO_NW*log(RHO_NW(-1)) + E_RHO_NW*std_E_RHO_NW ; % Price markup shock (THETAH_NW) log(THETAH_NW) = (1-LAMBDA_THETAH_NW)*log(THETAH_NW_SS) + LAMBDA_THETAH_NW*log(THETAH_NW(-1)) + E_THETAH_NW*std_E_THETAH_NW ; % Marginal efficiency of investment shock (Z_I_NW) log(Z_I_NW) = (1-LAMBDA_I_NW)*log(Z_I_NW_SS) + LAMBDA_I_NW*log(Z_I_NW(-1)) + E_I_NW*std_E_I_NW ; % Temporary labor augmenting technology shock (Z_L_NW) log(Z_L_NW) = (1-LAMBDA_L_NW)*log(Z_L_NW_SS) + LAMBDA_L_NW*log(Z_L_NW(-1)) + E_L_NW*std_E_L_NW ; % Consumption preference shock (Z_U_NW) log(Z_U_NW) = (1-LAMBDA_U_NW)*log(Z_U_NW_SS) + LAMBDA_U_NW*log(Z_U_NW(-1)) + E_U_NW*std_E_U_NW ; % Monetary policy shock (Z_RN3M_NW) log(Z_RN3M_NW) = (1-LAMBDA_RN3M_NW)*log(Z_RN3M_NW_SS) + LAMBDA_RN3M_NW*log(Z_RN3M_NW(-1)) + E_RN3M_NW*std_E_RN3M_NW; % Unit root processes Z/Z(-1) = DZT_NW_SS^(1-LAMBDA_DZT_NW)* (Z(-1)/Z(-2))^LAMBDA_DZT_NW*exp(std_E_DZT_NW*E_DZT_NW); UPSILON/UPSILON(-1) = DUT_NW_SS^(1-LAMBDA_DUT_NW)* 53 (UPSILON(-1)/UPSILON(-2))^LAMBDA_DUT_NW*exp(std_E_DUT_NW*E_DUT_NW); % Construct reported variables that can be asked for in IRFs, but % are not part of the model. Here you can use all MATLAB functions % that act on a double vector, and returns an output with the same % size as the input. reporting C_NW_LEVEL = exp(cumsum(log(DPQ_C_NW) - log(steady_state(DPQ_C_NW)))); I_NW_LEVEL = exp(cumsum(log(DPQ_I_NW) - log(steady_state(DPQ_I_NW)))); K_NW_LEVEL = exp(cumsum(log(DPQ_K_NW) - log(steady_state(DPQ_K_NW)))); Y_NW_LEVEL = exp(cumsum(log(DPQ_Y_NW) - log(steady_state(DPQ_Y_NW)))); The set of NB toolbox commands to transform and solve the model: % Read the non−stationary model modelNS = nb_dsge('nb_file','jpt_non_stationary.nb'); % Give it a name modelNS = set(modelNS,'name','Stationarize JPT automatically'); % Set the parameters param = struct(); param.ALPHA_NW = 0.167; param.BC_NW = 0.859; param.BETA_NW = 100/(0.134+100); param.BL_NW = 0; param.DELTA_NW = 0.025; param.DPQ_P_NW_SS = (0.702 + 100)/100; param.DUT_NW_SS = 1 + (0.597/100); param.DZT_NW_SS = 1 + (0.303 − (param.ALPHA_NW/... (1 − param.ALPHA_NW))*0.597)/100; param.LAMBDA_DUT_NW = 0.156; param.LAMBDA_DZT_NW = 0.286; param.LAMBDA_I_NW = 0.772; param.LAMBDA_L_NW = 0; param.LAMBDA_PSI_NW = 0.967; param.LAMBDA_RHO_NW = 0.590; param.LAMBDA_RN3M_NW = 0; param.LAMBDA_U_NW = 0; param.LAMBDA_THETAH_NW = 0.971; param.OMEGA_DPQ_Y_NW = 0.208; param.OMEGA_P_NW = 1.709; param.OMEGA_R_NW = 0.858; param.OMEGA_Y_NW = 0.051; param.PHI_PQ_NW = 0.2; param.PHI_I1_NW = 2.657; param.PHI_W_NW = 1.0080; param.PHI_U_NW = 5.434; param.PSI_NW_SS = 1.135/(1.135 − 1); param.RHO_NW_SS = 1; param.THETAH_NW_SS = 1.171/(1.171 − 1); param.ZETA_NW = 4.444; param.Z_I_NW_SS = 1; 54 % Z_U_NW shock process unitrootvars Z X model % 1) Final goods sector % Production function (A_NW) A_NW = Q_NW; % FOC (Q_NW) REAL_PQ_NW = 1; % 2) Intermediate goods sector % Intermidate production function (KBAR_NW) T_NW = (Z*Z_L_NW*L_NW)^(1-ALPHA_NW)*KBAR_NW^ALPHA_NW; % Optimality condition wrt utilized capital (I.e. demand function) (MC_NW) KBAR_NW = ALPHA_NW*(MC_NW/RK_NW)*T_NW; % Optimality condition wrt aggregated labor (I.e. demand function) (L_NW) L_NW = (1 - ALPHA_NW)*(MC_NW/REAL_W_NW)*T_NW; % Pricing (REAL_PQ_NW) Q_NW - THETAH_NW*Q_NW + MC_NW*THETAH_NW*Q_NW/REAL_PQ_NW - 100*PHI_PQ_NW*(DPQ_PQ_NW/DPQ_PQ_NW(-1) - 1)*DPQ_PQ_NW/DPQ_PQ_NW(-1)*Q_NW + DSA_NW*100*PHI_PQ_NW*(DPQ_PQ_NW(+1)/DPQ_PQ_NW - 1) *DPQ_PQ_NW(+1)^2/DPQ_PQ_NW*Q_NW(+1) = 0; % 3) Investment producer % Investment production function (Y_I_NW) I_NW = Y_I_NW; % First order condition for investment (REAL_PI_NW) REAL_PI_NW = 1; % 4) Capital producer % Capital production function (KNEW_NW) KNEW_NW = Z_I_NW*(1 - S_NW)*I_NW; % Investment adjustment cost function (S_NW) S_NW = (PHI_I1_NW/2)*(I_NW/I_NW(-1) - bgp(I_NW))^2; % Derivative of the investment adjustment cost function with respect to % first input (Multiplied by I_NW!) (S_PRIME1_NW) S_PRIME1_NW = PHI_I1_NW*(I_NW/I_NW(-1) - bgp(I_NW))*I_NW/I_NW(-1); % Derivative of the investment adjustment cost function with respect to % second input (Multiplied by I_NW!) (S_PRIME2_NW) S_PRIME2_NW = -PHI_I1_NW*(I_NW/I_NW(-1) - bgp(I_NW))*(I_NW/I_NW(-1))^2; 61 % Optimal capital investment (I_NW) 0 = DSA_NW*REAL_PK_NW(+1)*Z_I_NW(+1)*S_PRIME2_NW(+1) + REAL_PK_NW*Z_I_NW*(1 - S_NW - S_PRIME1_NW) - REAL_PI_NW; % 5) Households % Marignal utility of consumption (UPRIME_NW) UPRIME_NW = X*Z_U_NW*((C_NW - C_NW(-1)*BC_NW)/(1 - BC_NW/bgp(C_NW)))^(-1) ; % Marginal utility of labor labor (VPRIME_NW) VPRIME_NW = ((L_NW - BL_NW*L_NW(-1))/(1-BL_NW))^ZETA_NW ; % Stochastic discount factor (DSA_NW) DSA_NW = BETA_NW*RHO_NW(+1)*UPRIME_NW(+1)/(DPQ_P_NW(+1)*RHO_NW*UPRIME_NW); % Consumption euler equation (RN3M_NW) DSA_NW*RN3M_NW = 1; % Optimal wage setting (REAL_W_NW) REAL_W_NW = PSI_NW*MRS_NW/((PSI_NW-1)*(1-GAMMA_W_NW) + 1000*PHI_W_NW*DPQ_W_NW/DPQ_W_NW(-1)*(DPQ_W_NW/DPQ_W_NW(-1) - 1) - DSA_NW*DPQ_W_NW(+1)*L_NW(+1)/L_NW*1000*PHI_W_NW*DPQ_W_NW(+1) /DPQ_W_NW*(DPQ_W_NW(+1)/DPQ_W_NW - 1)); % Wage adjusment cost (GAMMA_W_NW) GAMMA_W_NW = 1000*PHI_W_NW/2*(DPQ_W_NW/DPQ_W_NW(-1) - 1)^2 ; % Marginal rate of substitution between consumption and leisure (MRS_NW) MRS_NW = VPRIME_NW/UPRIME_NW; % 6) Entrepreneurs % optimality with respect to capital (REAL_PK_NW) REAL_PK_NW = DSA_NW*DPQ_P_NW(+1)*( REAL_PK_NW(+1)*(1 - DELTA_NW) + RK_NW(+1)*U_NW(+1) - GAMMA_U_NW(+1)); % optimality with respect to utilization (RK_NW) RK_NW = GAMMAPRIME_U_NW; % Cost of utilizing the capital (GAMMA_U_NW) GAMMA_U_NW = steady_state(RK_NW)/PHI_U_NW*(exp(PHI_U_NW*(U_NW - 1)) - 1) ; % Marginal cost of utilizing capital (GAMMAPRIME_U_NW) GAMMAPRIME_U_NW = steady_state(RK_NW)*exp(PHI_U_NW*(U_NW - 1)) ; % Capital accumulation (K_NW) K_NW = (1-DELTA_NW)*K_NW(-1) + KNEW_NW; % Capital utilization (U_NW) KBAR_NW = U_NW*K_NW(-1); % 7) Market clearing % Final good market clearing (C_NW) A_NW = C_NW + Y_I_NW; 62 % Intermediate good market clearing (T_NW) T_NW = Q_NW; % 8) Definitions % Definition of natural output NAT_Y_NW = A_NW; % Consumption growth (DPQ_C_NW) DPQ_C_NW = C_NW/C_NW(-1); % Investment growth (DPQ_I_NW) DPQ_I_NW = I_NW/I_NW(-1); % Capital growth (DPQ_K_NW) DPQ_K_NW = K_NW/K_NW(-1); % Definition of natural output growth (DPQ_A_NW) DPQ_Y_NW = NAT_Y_NW/NAT_Y_NW(-1); % Definition of intermidate good inflation (DPQ_PQ_NW) DPQ_PQ_NW = DPQ_P_NW*REAL_PQ_NW/REAL_PQ_NW(-1); % Real investment inflation (DPQ_REAL_PI_NW) DPQ_REAL_PI_NW = REAL_PI_NW/REAL_PI_NW(-1); % Real wage inflation (DPQ_REAL_W_NW) DPQ_REAL_W_NW = REAL_W_NW/REAL_W_NW(-1); % Definition of wage inflation (DPQ_W_NW) DPQ_W_NW = REAL_W_NW/REAL_W_NW(-1)*DPQ_P_NW; % 9) Taylor rule (DPQ_P_NW) RN3M_NW = RN3M_NW(-1)^OMEGA_R_NW*(steady_state(RN3M_NW) *(DPQ_P_NW/steady_state(DPQ_P_NW))^OMEGA_P_NW *(NAT_Y_NW/steady_state(NAT_Y_NW))^OMEGA_Y_NW *(DPQ_Y_NW/steady_state(DPQ_Y_NW))^OMEGA_DPQ_Y_NW) ^(1-OMEGA_R_NW)*(Z_RN3M_NW/steady_state(Z_RN3M_NW)); [static] DPQ_P_NW = DPQ_P_NW_SS; % 10) Shock processes % Labor market competition shock (PSI_NW) log(PSI_NW) = (1-LAMBDA_PSI_NW)*log(PSI_NW_SS) + LAMBDA_PSI_NW*log(PSI_NW(-1)) + E_PSI_NW*std_E_PSI_NW ; % Discount factor shock (RHO_NW) log(RHO_NW) = (1-LAMBDA_RHO_NW)*log(RHO_NW_SS) + LAMBDA_RHO_NW*log(RHO_NW(-1)) + E_RHO_NW*std_E_RHO_NW ; % Price markup shock (THETAH_NW) log(THETAH_NW) = (1-LAMBDA_THETAH_NW)*log(THETAH_NW_SS) + LAMBDA_THETAH_NW*log(THETAH_NW(-1)) 63 + E_THETAH_NW*std_E_THETAH_NW ; % Marginal efficiency of investment shock (Z_I_NW) log(Z_I_NW) = (1-LAMBDA_I_NW)*log(Z_I_NW_SS) + LAMBDA_I_NW*log(Z_I_NW(-1)) + E_I_NW*std_E_I_NW ; % Temporary labor augmenting technology shock (Z_L_NW) log(Z_L_NW) = (1-LAMBDA_L_NW)*log(Z_L_NW_SS) + LAMBDA_L_NW*log(Z_L_NW(-1)) + E_L_NW*std_E_L_NW ; % Consumption preference shock (Z_U_NW) log(Z_U_NW) = (1-LAMBDA_U_NW)*log(Z_U_NW_SS) + LAMBDA_U_NW*log(Z_U_NW(-1)) + E_U_NW*std_E_U_NW ; % Monetary policy shock (Z_RN3M_NW) log(Z_RN3M_NW) = (1-LAMBDA_RN3M_NW)*log(Z_RN3M_NW_SS) + LAMBDA_RN3M_NW*log(Z_RN3M_NW(-1)) + E_RN3M_NW*std_E_RN3M_NW; % Unit root processes Z/Z(-1) = DZT_NW_SS^(1-LAMBDA_DZT_NW)*(Z(-1)/Z(-2))^LAMBDA_DZT_NW*exp(std_E_DZT_NW*E_DZT_NW); X/X(-1) = DXT_NW_SS^(1-LAMBDA_DXT_NW)*(X(-1)/X(-2))^LAMBDA_DXT_NW*exp(std_E_DXT_NW*E_DXT_NW); % Construct reported variables that can be asked for in IRFs, but % are not part of the model. Here you can use all MATLAB functions % that act on a double vector, and returns an output with the same % size as the input. reporting C_NW_LEVEL = exp(cumsum(log(DPQ_C_NW) - log(steady_state(DPQ_C_NW)))); I_NW_LEVEL = exp(cumsum(log(DPQ_I_NW) - log(steady_state(DPQ_I_NW)))); K_NW_LEVEL = exp(cumsum(log(DPQ_K_NW) - log(steady_state(DPQ_K_NW)))); Y_NW_LEVEL = exp(cumsum(log(DPQ_Y_NW) - log(steady_state(DPQ_Y_NW)))); With this le and the following commands to make the model stationary and compute impulse responses. % Read the non−stationary model modelNS = nb_dsge('nb_file','jpt_non_stationary_pref.nb'); % Give it a name modelNS = set(modelNS,'name','Stationarize JPT automatically'); % Set the parameters param = struct(); param.ALPHA_NW = 0.167; param.BC_NW = 0.859; param.BETA_NW = 100/(0.134+100); param.BL_NW = 0; param.DELTA_NW = 0.025; param.DPQ_P_NW_SS = (0.702 + 100)/100; param.DXT_NW_SS = 1 + (0.597/100); param.DZT_NW_SS = 1 + (0.303 − (param.ALPHA_NW/... (1 − param.ALPHA_NW))*0.597)/100; 64 param.LAMBDA_DXT_NW = 0.156; param.LAMBDA_DZT_NW = 0.286; param.LAMBDA_I_NW = 0.772; param.LAMBDA_L_NW = 0; param.LAMBDA_PSI_NW = 0.967; param.LAMBDA_RHO_NW = 0.590; param.LAMBDA_RN3M_NW = 0; param.LAMBDA_U_NW = 0; param.LAMBDA_THETAH_NW = 0.971; param.OMEGA_DPQ_Y_NW = 0.208; param.OMEGA_P_NW = 1.709; param.OMEGA_R_NW = 0.858; param.OMEGA_Y_NW = 0.051; param.PHI_PQ_NW = 0.2; param.PHI_I1_NW = 2.657; param.PHI_W_NW = 1.0080; param.PHI_U_NW = 5.434; param.PSI_NW_SS = 1.135/(1.135 − 1); param.RHO_NW_SS = 1; param.THETAH_NW_SS = 1.171/(1.171 − 1); param.ZETA_NW = 4.444; param.Z_I_NW_SS = 1; param.Z_L_NW_SS = 1; param.Z_U_NW_SS = 1; param.Z_RN3M_NW_SS = 1; param.std_E_DXT_NW = 0.630; param.std_E_DZT_NW = 0.933; param.std_E_I_NW = 5.103; param.std_E_L_NW = 0; param.std_E_PSI_NW = 0.310; param.std_E_RHO_NW = 0.036; param.std_E_RN3M_NW = 0.210; param.std_E_THETAH_NW = 0.219; param.std_E_U_NW = 0; modelNS = assignParameters(modelNS,param); % Solve for the balanced growth path (steps 1−3 of the algorithm) modelNS = solveBalancedGrowthPath(modelNS); % Stationarize the model (step 4 of the algorithm) modelNS = stationarize(modelNS) % Solve stationary steady state (numerically) modelNS = checkSteadyState(modelNS,... 'solver','fsolve',... 'steady_state_solve', true,... 'steady_state_default', @ones); % Return the steady state solution as a cell matrix ss = getSteadyState(modelNS) % Return the solution to the balanced growth path as a cell matrix bgp = getBalancedGrowthPath(modelNS) % Write the stationary model to file (optional) 65 writeModel2File(modelNS,'stationarized_pref.nb') % Obtain the state space representation of the solution % y(t) = A y(t−1) + B eps(t) modelNS = solve(modelNS); % compute and plot impulse responses [~,~,plotter] = irf([modelNS,modelS],'periods', 60,'plotSS', true,'shocks',... {'E_DXT_NW'},... 'variables',{'C_NW','DPQ_P_NW','DPQ_W_NW','I_NW','K_NW',... 'L_NW','NAT_Y_NW','RN3M_NW'},... 'settings',{'legBox','off','legFontSize',18,'subPlotSize',[3,3],... 'figureTitle',false,'lookUpMatrix','lookUpMatrixJPT',... 'legends',{'Stationarize JPT automatically',... 'Stationarize JPT manually',... '',... 'Steady state'}}); nb_graphInfoStructGUI(plotter) H Generating results from the JPT model In the main paper we have calculated the steady state, balanced growth path and moments. In this section we present the code on how you can calculate these tables using NB toolbox. %% Parametrization param = struct(); param.ALPHA_NW = 0.167; param.BC_NW = 0.859; param.BETA_NW = 100/(0.134+100); param.BL_NW = 0; param.DELTA_NW = 0.025; param.DPQ_P_NW_SS = (0.702 + 100)/100; param.DUT_NW_SS = 1 + (0.597/100); param.DZT_NW_SS = 1 + (0.303 − (param.ALPHA_NW/... (1 − param.ALPHA_NW))*0.597)/100; param.LAMBDA_DUT_NW = 0.156; param.LAMBDA_DZT_NW = 0.286; param.LAMBDA_I_NW = 0.772; param.LAMBDA_L_NW = 0; param.LAMBDA_PSI_NW = 0.967; param.LAMBDA_RHO_NW = 0.590; param.LAMBDA_RN3M_NW = 0; param.LAMBDA_U_NW = 0; param.LAMBDA_THETAH_NW = 0.971; param.OMEGA_DPQ_Y_NW = 0.208; param.OMEGA_P_NW = 1.709; param.OMEGA_R_NW = 0.858; param.OMEGA_Y_NW = 0.051; param.PHI_PQ_NW = 0.2; param.PHI_I1_NW = 2.657; param.PHI_W_NW = 1.0080; 66 param.PHI_U_NW = 5.434; param.PSI_NW_SS = 1.135/(1.135 − 1); param.RHO_NW_SS = 1; param.THETAH_NW_SS = 1.171/(1.171 − 1); param.ZETA_NW = 4.444; param.Z_I_NW_SS = 1; param.Z_L_NW_SS = 1; param.Z_U_NW_SS = 1; param.Z_RN3M_NW_SS = 1; param.std_E_DUT_NW = 0.630; param.std_E_DZT_NW = 0.933; param.std_E_I_NW = 5.103; param.std_E_L_NW = 0; param.std_E_PSI_NW = 0.310; param.std_E_RHO_NW = 0.036; param.std_E_RN3M_NW = 0.210; param.std_E_THETAH_NW = 0.219; param.std_E_U_NW = 0; %% Read the non−stationary model modelNS = nb_dsge('nb_file','jpt_non_stationary.nb'); modelNS = set(modelNS,'name','JPT with investment specific technology'); %% Assign parameters modelNS = assignParameters(modelNS,param); %% Solve for the balanced growth path modelNS = solveBalancedGrowthPath(modelNS); %% Stationarize the non−stationary model modelNS = stationarize(modelNS); %% Solve steady state numerically ssInit = struct(... 'GAMMA_W_NW',0,... 'PSI_NW',param.PSI_NW_SS,... 'S_NW',0,... 'S_PRIME1_NW',0,... 'S_PRIME2_NW',0,... 'GAMMA_U_NW',0,... 'THETAH_NW',param.THETAH_NW_SS); modelNS = checkSteadyState(modelNS,... 'solver','fsolve',... 'steady_state_solve', true,... 'steady_state_init', ssInit,... 'steady_state_default', @ones); %% Solve stationary model modelNS = solve(modelNS); 67 %% Parametrization paramP = struct(); paramP.ALPHA_NW = 0.167; paramP.BC_NW = 0.859; paramP.BETA_NW = 100/(0.134+100); paramP.BL_NW = 0; paramP.DELTA_NW = 0.025; paramP.DPQ_P_NW_SS = (0.702 + 100)/100; paramP.DXT_NW_SS = 1 + (0.597/100); paramP.DZT_NW_SS = 1 + (0.303 − (param.ALPHA_NW/... (1 − param.ALPHA_NW))*0.597)/100; paramP.LAMBDA_DXT_NW = 0.156; paramP.LAMBDA_DZT_NW = 0.286; paramP.LAMBDA_I_NW = 0.772; paramP.LAMBDA_L_NW = 0; paramP.LAMBDA_PSI_NW = 0.967; paramP.LAMBDA_RHO_NW = 0.590; paramP.LAMBDA_RN3M_NW = 0; paramP.LAMBDA_U_NW = 0; paramP.LAMBDA_THETAH_NW = 0.971; paramP.OMEGA_DPQ_Y_NW = 0.208; paramP.OMEGA_P_NW = 1.709; paramP.OMEGA_R_NW = 0.858; paramP.OMEGA_Y_NW = 0.051; paramP.PHI_PQ_NW = 0.2; paramP.PHI_I1_NW = 2.657; paramP.PHI_W_NW = 1.0080; paramP.PHI_U_NW = 5.434; paramP.PSI_NW_SS = 1.135/(1.135 − 1); paramP.RHO_NW_SS = 1; paramP.THETAH_NW_SS = 1.171/(1.171 − 1); paramP.ZETA_NW = 4.444; paramP.Z_I_NW_SS = 1; paramP.Z_L_NW_SS = 1; paramP.Z_U_NW_SS = 1; paramP.Z_RN3M_NW_SS = 1; paramP.std_E_DXT_NW = 0.630; paramP.std_E_DZT_NW = 0.933; paramP.std_E_I_NW = 5.103; paramP.std_E_L_NW = 0; paramP.std_E_PSI_NW = 0.310; paramP.std_E_RHO_NW = 0.036; paramP.std_E_RN3M_NW = 0.210; paramP.std_E_THETAH_NW = 0.219; paramP.std_E_U_NW = 0; %% Read the non−stationary model modelNSPref = nb_dsge('nb_file','jpt_non_stationary_pref.nb'); modelNSPref = set(modelNSPref,'name','JPT with unit root in preferences'); %% Assign parameters modelNSPref = assignParameters(modelNSPref,paramP); 68 %% Solve for the balanced growth path modelNSPref = solveBalancedGrowthPath(modelNSPref); %% Stationarize the non−stationary model modelNSPref = stationarize(modelNSPref); %% Solve steady state numerically ssInit = struct(... 'GAMMA_W_NW',0,... 'PSI_NW',paramPref.PSI_NW_SS,... 'S_NW',0,... 'S_PRIME1_NW',0,... 'S_PRIME2_NW',0,... 'GAMMA_U_NW',0,... 'THETAH_NW',paramPref.THETAH_NW_SS); modelNSPref = checkSteadyState(modelNSPref,... 'solver','fsolve',... 'steady_state_solve', true,... 'steady_state_init', ssInit,... 'steady_state_default', @ones); %% Solve stationary model modelNSPref = solve(modelNSPref); %% Compare balanced growth paths ssGrowth = getBalancedGrowthPath([modelNS,modelNSPref],... {'C_NW','I_NW','K_NW','L_NW','NAT_Y_NW','REAL_W_NW',... 'REAL_PI_NW'},'headers') %% Compare steady state results ss = getSteadyState([modelNS,modelNSPref],... {'C_NW','DPQ_P_NW','DPQ_W_NW','I_NW','K_NW',... 'L_NW','NAT_Y_NW','RN3M_NW'},... 'headers') %% Theoretical moments vars = {'C_NW','DPQ_P_NW','DPQ_W_NW','I_NW','K_NW',... 'L_NW','NAT_Y_NW','RN3M_NW'}; [~,CNS] = theoreticalMoments(modelNS,'vars',vars,... 'type','covariance'); VNS = diag(CNS)'; VNS = rename(VNS,'variable','diag',getName(modelNS)); [~,CNSPref] = theoreticalMoments(modelNSPref,'vars',vars,... 'type','covariance'); VNSPref = diag(CNSPref)'; 69 VNSPref = rename(VNSPref,'variable','diag',getName(modelNSPref)); I The automatically made stationary JPT model le with unit root in preferences endogenous Z_U_NW Z_RN3M_NW Z_L_NW Z_I_NW Y_I_NW VPRIME_NW U_NW UPRIME_NW T_NW THETAH_NW S_PRIME2_NW S_PRIME1_NW S_NW RN3M_NW RK_NW RHO_NW REAL_W_NW REAL_PQ_NW REAL_PK_NW REAL_PI_NW Q_NW PSI_NW NAT_Y_NW MRS_NW MC_NW L_NW K_NW KNEW_NW KBAR_NW I_NW GAMMA_W_NW GAMMA_U_NW GAMMAPRIME_U_NW D_Z_Z D_Z_Y_I_NW D_Z_X D_Z_VPRIME_NW D_Z_UPRIME_NW D_Z_T_NW D_Z_REAL_W_NW D_Z_Q_NW D_Z_NAT_Y_NW D_Z_MRS_NW D_Z_L_NW D_Z_K_NW D_Z_KNEW_NW D_Z_KBAR_NW D_Z_I_NW D_Z_C_NW D_Z_A_NW DSA_NW DPQ_Y_NW DPQ_W_NW DPQ_REAL_W_NW DPQ_REAL_PI_NW DPQ_P_NW DPQ_PQ_NW DPQ_I_NW DPQ_C_NW C_NW A_NW exogenous E_DXT_NW E_DZT_NW E_I_NW E_L_NW E_PSI_NW E_RHO_NW E_RN3M_NW E_THETAH_NW E_U_NW parameters std_E_U_NW std_E_THETAH_NW std_E_RN3M_NW std_E_RHO_NW std_E_PSI_NW std_E_L_NW std_E_I_NW std_E_DZT_NW std_E_DXT_NW Z_U_NW_SS Z_RN3M_NW_SS Z_L_NW_SS Z_I_NW_SS ZETA_NW THETAH_NW_SS RHO_NW_SS PSI_NW_SS PHI_W_NW PHI_U_NW PHI_PQ_NW PHI_I1_NW OMEGA_Y_NW OMEGA_R_NW OMEGA_P_NW OMEGA_DPQ_Y_NW LAMBDA_U_NW LAMBDA_THETAH_NW LAMBDA_RN3M_NW LAMBDA_RHO_NW LAMBDA_PSI_NW LAMBDA_L_NW LAMBDA_I_NW LAMBDA_DZT_NW LAMBDA_DXT_NW DZT_NW_SS DXT_NW_SS DPQ_P_NW_SS DELTA_NW BL_NW BETA_NW BC_NW ALPHA_NW model A_NW-Q_NW; REAL_PQ_NW-1; T_NW = (((1*Z_L_NW)*L_NW)^(1-ALPHA_NW))*(KBAR_NW^ALPHA_NW); KBAR_NW = (ALPHA_NW*(MC_NW/RK_NW))*T_NW; L_NW = ((1-ALPHA_NW)*(MC_NW/REAL_W_NW))*T_NW; Q_NW-(THETAH_NW*Q_NW)+((MC_NW*THETAH_NW)*Q_NW)/REAL_PQ_NW- (((((100*PHI_PQ_NW)*(DPQ_PQ_NW/DPQ_PQ_NW(-1)-1))*DPQ_PQ_NW)/DPQ_PQ_NW(-1))*Q_NW)+ (((((DSA_NW*100)*PHI_PQ_NW)*(DPQ_PQ_NW(+1)/DPQ_PQ_NW-1))* (DPQ_PQ_NW(+1)^2))/DPQ_PQ_NW)*(Q_NW(+1)*D_Z_Q_NW(+1))-0; I_NW-Y_I_NW; REAL_PI_NW-1; KNEW_NW = (Z_I_NW*(1-S_NW))*I_NW; S_NW = (PHI_I1_NW/2)*((I_NW/(I_NW(-1)*D_Z_I_NW^-1)- steady_state(D_Z_I_NW))^2); S_PRIME1_NW = ((PHI_I1_NW*(I_NW/(I_NW(-1)*D_Z_I_NW^-1)- steady_state(D_Z_I_NW)))*I_NW)/(I_NW(-1)*D_Z_I_NW^-1); S_PRIME2_NW = ((-PHI_I1_NW)*(I_NW/(I_NW(-1)*D_Z_I_NW^-1)- steady_state(D_Z_I_NW)))*((I_NW/(I_NW(-1)*D_Z_I_NW^-1))^2); 0 = ((DSA_NW*REAL_PK_NW(+1))*Z_I_NW(+1))*S_PRIME2_NW(+1)+ 70