MATLAB programs by Peter Ireland

September 30, 2004

Contents 1 Introduction . . . . . . . . . . 2 Estimation . . . . . . . . . . . 2.1 est.m . . . . . . . . . . . 2.2 estse.m . . . . . . . . . . 2.3 estseq.m . . . . . . . . . 2.4 llfn.m . . . . . . . . . . 2.5 llfnse.m . . . . . . . . . 3 Solution and Simulation . . . 3.1 solv.m . . . . . . . . . . 3.2 check.m . . . . . . . . . 3.3 ksmooth.m . . . . . . . . 4 Impulse Responses & Variance 4.1 imp.m . . . . . . . . . . 4.2 vardec.m . . . . . . . . . 4.3 vardecfn.m . . . . . . . . 5 Stability Tests . . . . . . . . . 5.1 stabtest.m . . . . . . . . 6 Forecasting . . . . . . . . . . 6.1 fork.m . . . . . . . . . . 6.2 forkfn.m . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Decompostion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . .

1 3 3 9 13 19 23 27 27 29 31 35 35 37 39 45 45 52 52 60

1. Introduction Peter Ireland has made MATLAB programs (i.e. both scripts and functions) available so the econometric work in his paper ”A Method for Taking Models to the Data,” in the March 2004 issue of the Journal of Economic Dynamics and Control can be replicated. These programs can be found at: http://www2.bc.edu/~irelandp/#programs. A brief taxonomy of the programs according to the tasks performed include:

Estimation • est.m: Maximizes (minimizes negative) log-likelihood function for Hansen’s RBC model to obtain the parameter estimates where some have been transformed to comply with theoretical restrictions. • estse.m: Computes standard errors for the transformed parameter estimates. • estseq.m: Recursive estimation of est.m with variable end dates. • llfn.m: MATLAB function called by est.m and estseq.m. Uses the Kalman filter to evaluate the log-likelihood function based on the transformed parameters. • llfnse.m: MATLAB function called by estse.m. Uses the Kalman filter to evaluate the log-likelihood function based on the untransformed parameters.

Solution & Simulation • solv.m: Solves Hansen’s RBC model. • check.m: Checks the solution found by the script solv.m. • ksmooth.m: Provides smoothed estimates of the shocks in Hansen’s model (i.e. technology plus measurement errors). Also calculates the correlation between the technology shock and the measurement errors.

Impulse Responses & Variance Decomposition • imp.m: Computes impulse responses for Hansen’s model using the solutions provided by solv.m.

1

• vardec.m: Uses the parameter estimates and standard errors provided by estse.m to calculate variance decompositions and standard errors for Hansen’s model. • vardecfn.m: MATLAB function which does the same as vardec.m (less the standard errors) for a given set of parameter estimates.

Stability Tests stabtest.m: Performs parameter stability tests for the Hansen model (break point occurs at 1973:1).

Forecasting • fork.m: Generates k-step-ahead forecasts from four competing models: Hansen’s model with vector AR(1) and scalar AR(1) errors plus a UVAR(1) and a UVAR(2) model. • forkfn.m: MATLAB function which uses the Kalman filter to generate k-stepahead forecasts for the non-predetermined variables in Hansen’s model.

Extra MATLAB Programs • The MATLAB programs briefly described above and fully listed below are relevant for the estimation and analysis of Ireland’s full model (i.e. with vector AR(1) measurement errors). His model can also be estimated with scalar AR(1) errors (i.e. by assuming that the covariance matrix of his measurement equations, D and the covariance matrix of his measurement equations errors, V2 are diagonal). His MATLAB scripts estdse.m and estseqd.m perform the same tasks as the above described programs estse.m and estseq.m except V2 and D are diagonal. This also applies to the MATLAB functions llfnd.m and llfndse.m which correspond to llfn.m and llfnse.m. The program stabtst2.m performs the same stability tests as stabtest.m except the break point is 1980:1. Finally the functions llfn2.m & llfnse2.m and llfn3.m & llfnse3.m correspond to llfn.m & llfnse.m except the time trend is adjusted to accommodate the post 1973 and post 1981 sub-samples respectively.

2

2. Estimation 2.1. est.m % est.m: Maximizes (minimizes negative) log likelihood function % for the real business cycle model with indivisible labor. % % When maximizing the log likelihood function, the parameters % are transformed so that they satisfy theoretical % restrictions. The log likelihood function with % transformed parameters is contained in llfn.m. % % T H IS P R O G R A M WA S W R IT T E N FO R M AT L A B B Y P E T E R N . IR E L A N D % B O S T O N C O L L E G E , D E PA RT M E N T O F E C O N O M IC S 1 4 0 C O M M O N W E A LT H AV E N U E , % C H E S T N U T H IL L , M A 0 2 4 67 , irela n d p @ b c.e d u % % F IN A N C IA L S U P P O RT F R O M T H E N AT IO N A L S C IE N C E FO U N D AT IO N U N D E R G R A N T N O S . % S E S -9 9 8 57 6 3 A N D S E S -0 2 1 3 4 6 1 IS G R AT E F U L LY A C K N O W L E D G E D . % % C O P Y R IG H T (c) 20 0 3 B Y P E T E R N . IR E L A N D . R E % D IS T R IB U T IO N IS P E R M IT T E D FO R % E D U C AT IO N A L A N D R E S E A R C H P U R P O S E S , S O L O N G A S N O C H A N G E S A R E M A D E . A L L % C O P IE S M U S T B E P R O V ID E D F R E E O F C H A R G E A N D M U S T IN C L U D E T H IS C O P Y R IG H T N O T IC E .

% load data and choose sample period global yt ct ht load ych.dat; % full sample yt = ych(:,1); ct = ych(:,2); ht = ych(:,3); % pre-1973 subsample % yt = ych(1:100,1); % ct = ych(1:100,2); % ht = ych(1:100,3); % post-1973 subsample % yt = ych(101:218,1); % ct = ych(101:218,2); 3

% ht = ych(101:218,3); % pre-1980 subsample % yt = ych(1:128,1); % ct = ych(1:128,2); % ht = ych(1:128,3); % post-1980 subsample % yt = ych(129:218,1); % ct = ych(129:218,2); % ht = ych(129:218,3); % set starting values for full sample bettr = sqrt(0.99/(1-0.99)); gamtr = 0.0045; thettr = sqrt(0.20/(1-0.20)); etatr = 0.0051; delttr = sqrt(0.025/(1-0.025)); atr = 6; rhotr = 0.9975; sigtr = 0.0055; dyytr = 1.4187; dyctr = 0.2251; dyhtr = -0.4441; dcytr = 0.0935; dcctr = 1.0236; dchtr = -0.0908; dhytr = 0.7775; dhctr = 0.3706; dhhtr = 0.2398; vyytr = 0.0072; vcytr = 0.0040; vcctr = 0.0057; vhytr = 0.0015; vhctr = 0.0010; vhhtr = 0.0000; % set starting values for pre-1973 sample % bettr = sqrt(0.99/(1-0.99)); 4

% gamtr = 0.0045; % thettr = sqrt(0.20/(1-0.20)); % etatr = 0.0051; % delttr = sqrt(0.025/(1-0.025)); % atr = 6; % rhotr = 0.9964; % sigtr = 0.0054; % dyytr = 1.1883; % dyctr = 0.3612; % dyhtr = -0.4075; % dcytr = 0.0875; % dcctr = 1.0024; % dchtr = -0.0738; % dhytr = 0.6527; % dhctr = 0.3161; % dhhtr = 0.3644; % vyytr = 0.0090; % vcytr = 0.0061; % vcctr = -0.0064; % vhytr = 0.0029; % vhctr = -0.0011; % vhhtr = -0.0031; % set starting values for post-1973 sample % bettr = sqrt(0.99/(1-0.99)); % gamtr = 0.0045; % thettr = sqrt(0.20/(1-0.20)); % etatr = 0.0051; % delttr = sqrt(0.025/(1-0.025)); % atr = 6; % rhotr = 0.9991; % sigtr = 0.0048; % dyytr = 1.2245; % dyctr = 0.7933; % dyhtr = -0.4947; % dcytr = 0.0285; 5

% dcctr = 1.0665; % dchtr = -0.0560; % dhytr = 0.3008; % dhctr = 0.9692; % dhhtr = 0.3935; % vyytr = 0.0057; % vcytr = 0.0034; % vcctr = -0.0041; % vhytr = 0.0012; % vhctr = -0.0015; % vhhtr = 0.0000; % set starting values for pre-1980 sample % bettr = sqrt(0.99/(1-0.99)); % gamtr = 0.0045; % thettr = sqrt(0.20/(1-0.20)); % etatr = 0.0051; % delttr = sqrt(0.025/(1-0.025)); % atr = 6; % rhotr = 0.9959; % sigtr = 0.0052; % dyytr = 1.2635; % dyctr = 0.2386; % dyhtr = -0.4115; % dcytr = 0.0710; % dcctr = 1.0475; % dchtr = -0.1075; % dhytr = 0.6708; % dhctr = 0.2053; % dhhtr = 0.3779; % vyytr = -0.0087; % vcytr = -0.0056; % vcctr = 0.0065; % vhytr = -0.0030; % vhctr = 0.0006; % vhhtr = -0.0034; 6

% set starting values for post-1980 sample % bettr = sqrt(0.99/(1-0.99)); % gamtr = 0.0045; % thettr = sqrt(0.20/(1-0.20)); % etatr = 0.0051; % delttr = sqrt(0.025/(1-0.025)); % atr = 6; % rhotr = 0.9985; % sigtr = 0.0041; % dyytr = 1.2801; % dyctr = 0.4877; % dyhtr = -0.4704; % dcytr = 0.0533; % dcctr = 1.0441; % dchtr = -0.0714; % dhytr = 0.3816; % dhctr = 0.6283; % dhhtr = 0.3914; % vyytr = -0.0062; % vcytr = -0.0026; % vcctr = -0.0042; % vhytr = -0.0016; % vhctr = -0.0007; % vhhtr = 0.0000; % maximize likelihood bigtheto = [ gamtr thettr etatr atr rhotr sigtr ... dyytr dyctr dyhtr ... dcytr dcctr dchtr ... dhytr dhctr dhhtr ... vyytr ... vcytr vcctr ... vhytr vhctr vhhtr ]’; options(1) = 1; options(14) = 10000; thetstar = fminu(’llfn’,bigtheto,options); 7

% thetstar = fminu(’llfn2’,bigtheto,options); % thetstar = fminu(’llfn3’,bigtheto,options);

8

2.2. estse.m % estse.m: Maximizes (minimizes negative) log likelihood function % for the real business cycle model with indivisible labor % and computes standard errors for estimated parameters. % % When maximizing the log likelihood function, the parameters % are transformed so that they satisfy theoretical % restrictions. The log likelihood function with % transformed parameters is contained in llfn.m. % % When calculating standard errors, the original % parameters are used. The log likelihood function without % transformed parameters is contained in llfnse.m. % % The untransformed parameter estimates are reported as % the vector tstar. The standard errors are reported as % the vector sevec. % % T H IS P R O G R A M WA S W R IT T E N FO R M AT L A B B Y P E T E R N . IR E L A N D % B O S T O N C O L L E G E , D E PA RT M E N T O F E C O N O M IC S 1 4 0 C O M M O N W E A LT H AV E N U E , % C H E S T N U T H IL L , M A 0 2 4 67 , irela n d p @ b c.e d u % % F IN A N C IA L S U P P O RT F R O M T H E N AT IO N A L S C IE N C E FO U N D AT IO N U N D E R G R A N T N O S . % S E S -9 9 8 57 6 3 A N D S E S -0 2 1 3 4 6 1 IS G R AT E F U L LY A C K N O W L E D G E D . % % C O P Y R IG H T (c) 20 0 3 B Y P E T E R N . IR E L A N D . R E % D IS T R IB U T IO N IS P E R M IT T E D FO R % E D U C AT IO N A L A N D R E S E A R C H P U R P O S E S , S O L O N G A S N O C H A N G E S A R E M A D E . A L L % C O P IE S M U S T B E P R O V ID E D F R E E O F C H A R G E A N D M U S T IN C L U D E T H IS C O P Y R IG H T N O T IC E .

% load data global yt ct ht scalinv load ych.dat; yt = ych(:,1); ct = ych(:,2); ht = ych(:,3); % set starting values 9

bettr = sqrt(0.99/(1-0.99)); gamtr = 0.0045; thettr = sqrt(0.20/(1-0.20)); etatr = 0.0051; delttr = sqrt(0.025/(1-0.025)); atr = 6; rhotr = 0.9975; sigtr = 0.0055; dyytr = 1.4187; dyctr = 0.2251; dyhtr = -0.4441; dcytr = 0.0935; dcctr = 1.0236; dchtr = -0.0908; dhytr = 0.7775; dhctr = 0.3706; dhhtr = 0.2398; vyytr = 0.0072; vcytr = 0.0040; vcctr = 0.0057; vhytr = 0.0015; vhctr = 0.0010; vhhtr = 0.0000; % maximize likelihood bigtheto = [ gamtr thettr etatr atr rhotr sigtr ... dyytr dyctr dyhtr ... dcytr dcctr dchtr ... dhytr dhctr dhhtr ... vyytr ... vcytr vcctr ... vhytr vhctr vhhtr ]’; options(1) = 1; options(14) = 10000; thetstar = fminu(’llfn’,bigtheto,options); % find standard errors 10

thetstar = real(thetstar); beta = bettr^2/(1+bettr^2); gamma = abs(thetstar(1)); theta = thetstar(2)^2/(1+thetstar(2)^2); eta = 1 + abs(thetstar(3)); delta = delttr^2/(1+delttr^2); a = abs(thetstar(4)); rho = thetstar(5); sig = abs(thetstar(6)); dyy = thetstar(7); dyc = thetstar(8); dyh = thetstar(9); dcy = thetstar(10); dcc = thetstar(11); dch = thetstar(12); dhy = thetstar(13); dhc = thetstar(14); dhh = thetstar(15); cholv = [ thetstar(16) 0 0 ; ... thetstar(17:18)’ 0 ; ... thetstar(19:21)’ ]; bigv = cholv*cholv’; vyy = sqrt(bigv(1,1)); vcc = sqrt(bigv(2,2)); vhh = sqrt(bigv(3,3)); vyc = bigv(1,2); vyh = bigv(1,3); vch = bigv(2,3); tstar = [ gamma theta eta a rho sig ... dyy dyc dyh dcy dcc dch dhy dhc dhh ... vyy vcc vhh vyc vyh vch ]’; scalvec = ones(21,1); scalinv = inv(diag(scalvec)); tstars = diag(scalvec)*tstar; fstar = llfnse(tstars); 11

eee = 1e-6; epsmat = eee*eye(21); hessvec = zeros(21,1); for i = 1:21 hessvec(i) = llfnse(tstars+epsmat(:,i)); end hessmat = zeros(21,21); for i = 1:21 for j = 1:21 hessmat(i,j) = (llfnse(tstars+epsmat(:,i)+epsmat(:,j)) ... -hessvec(i)-hessvec(j)+fstar)/eee^2; end end bighx = scalinv*inv(hessmat)*scalinv’; sevec = sqrt(diag(bighx));

12

2.3. estseq.m % estseq.m: Estimates the parameters of the real business cycle % model with indivisible labor for a sequence of sample periods % of increasing length. Starts with the sample period % 1948:1-1984:4 and ends with the sample period % 1948:1-2002:2. Thus, estimates the parameters for a % total of 71 samples. % % Returns: % thetmat = columns contain transformed estimates for each % sample % tstarmat = columns contain untransformed estimates for % each sample % % T H IS P R O G R A M WA S W R IT T E N FO R M AT L A B B Y P E T E R N . IR E L A N D % B O S T O N C O L L E G E , D E PA RT M E N T O F E C O N O M IC S 1 4 0 C O M M O N W E A LT H AV E N U E , % C H E S T N U T H IL L , M A 0 2 4 67 , irela n d p @ b c.e d u % % F IN A N C IA L S U P P O RT F R O M T H E N AT IO N A L S C IE N C E FO U N D AT IO N U N D E R G R A N T N O S . % S E S -9 9 8 57 6 3 A N D S E S -0 2 1 3 4 6 1 IS G R AT E F U L LY A C K N O W L E D G E D . % % C O P Y R IG H T (c) 20 0 3 B Y P E T E R N . IR E L A N D . R E % D IS T R IB U T IO N IS P E R M IT T E D FO R % E D U C AT IO N A L A N D R E S E A R C H P U R P O S E S , S O L O N G A S N O C H A N G E S A R E M A D E . A L L % C O P IE S M U S T B E P R O V ID E D F R E E O F C H A R G E A N D M U S T IN C L U D E T H IS C O P Y R IG H T N O T IC E .

% load data and set up output matrix global yt ct ht load ych.dat; bigt = length(ych(:,1)); thetmat = zeros(21,71); tstarmat = zeros(21,71); % set starting values gamtr = 0.0045; thettr = sqrt(0.20/(1-0.20)); etatr = 0.0051; atr = 6; 13

rhotr = 0.9975; sigtr = 0.0055; dyytr = 1.4187; dyctr = 0.2251; dyhtr = -0.4441; dcytr = 0.0935; dcctr = 1.0236; dchtr = -0.0908; dhytr = 0.7775; dhctr = 0.3706; dhhtr = 0.2398; vyytr = 0.0072; vcytr = 0.0040; vcctr = 0.0057; vhytr = 0.0015; vhctr = 0.0010; vhhtr = 0.0000; bigthet1 = [ gamtr thettr etatr atr rhotr sigtr ... dyytr dyctr dyhtr ... dcytr dcctr dchtr ... dhytr dhctr dhhtr ... vyytr ... vcytr vcctr ... vhytr vhctr vhhtr ]’; gamtr = 0.0045; thettr = 0.5441; etatr = 0.0050; atr = 5.2174; rhotr = 0.9986; sigtr = 0.0057; dyytr = 1.3897; dyctr = 0.3617; dyhtr = -0.5032; dcytr = 0.1471; dcctr = 0.9692; 14

dchtr = -0.1138; dhytr = 0.7510; dhctr = 0.4285; dhhtr = 0.2050; vyytr = 0.0071; vcytr = 0.0044; vcctr = 0.0056; vhytr = 0.0014; vhctr = 0.0012; vhhtr = 0.0000; bigthet2 = [ gamtr thettr etatr atr rhotr sigtr ... dyytr dyctr dyhtr ... dcytr dcctr dchtr ... dhytr dhctr dhhtr ... vyytr ... vcytr vcctr ... vhytr vhctr vhhtr ]’; gamtr = 0.0046; thettr = 0.5367; etatr = 0.0048; atr = 5.3634; rhotr = 0.9987; sigtr = 0.0058; dyytr = 1.3672; dyctr = 0.3997; dyhtr = -0.5324; dcytr = 0.1404; dcctr = 0.9805; dchtr = -0.1230; dhytr = 0.7083; dhctr = 0.4577; dhhtr = 0.2097; vyytr = 0.0070; vcytr = 0.0046; vcctr = 0.0056; 15

vhytr = 0.0012; vhctr = 0.0014; vhhtr = 0.0000; bigthet3 = [ gamtr thettr etatr atr rhotr sigtr ... dyytr dyctr dyhtr ... dcytr dcctr dchtr ... dhytr dhctr dhhtr ... vyytr ... vcytr vcctr ... vhytr vhctr vhhtr ]’; % estimate parameters for each sample period for t = bigt:-1:148 yt = ych(1:t,1); ct = ych(1:t,2); ht = ych(1:t,3); options(1) = 1; options(14) = 100000; if t == bigt thetstar = fminu(’llfn’,bigthet1,options); elseif t == 217 thetstar = fminu(’llfn’,bigthet1,options); elseif t == 210 thetstar = fminu(’llfn’,bigthet1,options); elseif t == 207 thetstar = fminu(’llfn’,bigthet2,options); elseif t == 206 thetstar = fminu(’llfn’,bigthet2,options); elseif t == 205 thetstar = fminu(’llfn’,bigthet2,options); elseif t == 201 thetstar = fminu(’llfn’,bigthet2,options); elseif t == 199 thetstar = fminu(’llfn’,bigthet2,options); elseif t == 195 thetstar = fminu(’llfn’,bigthet2,options); 16

elseif t == 192 thetstar = fminu(’llfn’,bigthet3,options); elseif t == 191 thetstar = fminu(’llfn’,bigthet3,options); elseif t == 188 thetstar = fminu(’llfn’,bigthet3,options); elseif t == 186 thetstar = fminu(’llfn’,bigthet3,options); elseif t == 183 thetstar = fminu(’llfn’,bigthet3,options); elseif t == 179 thetstar = fminu(’llfn’,bigthet3,options); elseif t == 172 thetstar = fminu(’llfn’,bigthet3,options); elseif t == 158 thetstar = fminu(’llfn’,bigthet3,options); elseif t == 157 thetstar = fminu(’llfn’,bigthet3,options); elseif t == 152 thetstar = fminu(’llfn’,bigthet3,options); else thetstar = fminu(’llfn’,thetstar,options); end thetstar = real(thetstar); gamma = abs(thetstar(1)); theta = thetstar(2)^2/(1+thetstar(2)^2); eta = 1 + abs(thetstar(3)); a = abs(thetstar(4)); rho = thetstar(5); sig = abs(thetstar(6)); dyy = thetstar(7); dyc = thetstar(8); dyh = thetstar(9); dcy = thetstar(10); dcc = thetstar(11); 17

dch = thetstar(12); dhy = thetstar(13); dhc = thetstar(14); dhh = thetstar(15); cholv = [ thetstar(16) 0 0 ; ... thetstar(17:18)’ 0 ; ... thetstar(19:21)’ ]; bigv = cholv*cholv’; vyy = sqrt(bigv(1,1)); vcc = sqrt(bigv(2,2)); vhh = sqrt(bigv(3,3)); vyc = bigv(1,2); vyh = bigv(1,3); vch = bigv(2,3); tstar = [ gamma theta eta a rho sig ... dyy dyc dyh dcy dcc dch dhy dhc dhh ... vyy vcc vhh vyc vyh vch ]’; thetmat(:,t-147) = thetstar; tstarmat(:,t-147) = tstar; end % save results save thetmat thetmat save tstarmat tstarmat

18

2.4. llfn.m function llfn = llfn(bigthet); % Uses the Kalman filter to evaluate the negative log likelihood % function for the real business cycle model with indivisible % labor. Some parameters are transformed so that they satisty % theoretical restrictions. % % T H IS P R O G R A M WA S W R IT T E N FO R M AT L A B B Y P E T E R N . IR E L A N D % B O S T O N C O L L E G E , D E PA RT M E N T O F E C O N O M IC S 1 4 0 C O M M O N W E A LT H AV E N U E , % C H E S T N U T H IL L , M A 0 2 4 67 , irela n d p @ b c.e d u % % F IN A N C IA L S U P P O RT F R O M T H E N AT IO N A L S C IE N C E FO U N D AT IO N U N D E R G R A N T N O S . % S E S -9 9 8 57 6 3 A N D S E S -0 2 1 3 4 6 1 IS G R AT E F U L LY A C K N O W L E D G E D . % % C O P Y R IG H T (c) 20 0 3 B Y P E T E R N . IR E L A N D . R E % D IS T R IB U T IO N IS P E R M IT T E D FO R % E D U C AT IO N A L A N D R E S E A R C H P U R P O S E S , S O L O N G A S N O C H A N G E S A R E M A D E . A L L % C O P IE S M U S T B E P R O V ID E D F R E E O F C H A R G E A N D M U S T IN C L U D E T H IS C O P Y R IG H T N O T IC E .

% define variables and parameters global yt ct ht bigt = length(yt); bigthet = real(bigthet); bettr = sqrt(0.99/(1-0.99)); gamtr = bigthet(1); thettr = bigthet(2); etatr = bigthet(3); delttr = sqrt(0.025/(1-0.025)); atr = bigthet(4); rhotr = bigthet(5); sigtr = bigthet(6); dyytr = bigthet(7); dyctr = bigthet(8); dyhtr = bigthet(9); dcytr = bigthet(10); dcctr = bigthet(11); dchtr = bigthet(12); 19

dhytr = bigthet(13); dhctr = bigthet(14); dhhtr = bigthet(15); vyytr = bigthet(16); vcytr = bigthet(17); vcctr = bigthet(18); vhytr = bigthet(19); vhctr = bigthet(20); vhhtr = bigthet(21); % untransform parameters beta = bettr^2/(1+bettr^2); gamma = abs(gamtr); theta = thettr^2/(1+thettr^2); eta = 1 + abs(etatr); delta = delttr^2/(1+delttr^2); a = abs(atr); rho = rhotr; sig = abs(sigtr); % compute steady-state values kappa = eta/beta - 1 + delta; lambda = eta - 1 + delta; hss = ((1-theta)/gamma)/(1-theta*lambda/kappa); yss = (a^(1/(1-theta)))*((theta/kappa)^(theta/(1-theta)))*hss; kss = (theta/kappa)*yss; iss = (theta*lambda/kappa)*yss; css = yss - iss; % compute K coefficients bigk11 = (eta-beta*(1-theta)*(1-delta))/(beta*eta*theta); bigk12 = (beta*eta*theta^2-eta+beta*(1-theta^2)*(1-delta)) ... /(beta*eta*theta^2); bigk22 = eta*theta/(eta-beta*(1-theta)*(1-delta)); % compute L coefficients bigl1 = (eta-beta*(1-delta))/(beta*eta*theta^2); bigl2 = rho*(eta-beta*(1-delta))/(eta-beta*(1-theta)*(1-delta)); % form S matrices 20

bigs1 = (bigk22-bigk11)/bigk12; bigs2 = ((bigk22-bigk11)*bigl1-bigk12*bigl2)/(bigk12*(bigk11-rho)); bigs3 = bigk22; bigs4 = bigk12*bigl2/(bigk22-bigk11) ... + (bigk22-rho)*((bigk22-bigk11)*bigl1-bigk12*bigl2) ... / ((bigk22-bigk11)*(bigk11-rho)); bigs5 = [ 1 - ((1-theta)/theta)*bigs1 ; ... bigs1 + (kappa/(theta^2*lambda))*(theta-bigs1) ; ... 1 - (1/theta)*bigs1 ]; bigs6 = [ 1/theta - ((1-theta)/theta)*bigs2 ; ... bigs2 + (kappa/(theta^2*lambda))*(1-bigs2) ; ... 1/theta - (1/theta)*bigs2 ]; % form matrices PI, W, and U bigpi = [ bigs3 bigs4 ; 0 rho ]; bigw = [ 0 ; 1 ]; bigu = [ bigs5 bigs6 ; bigs1 bigs2 ]; % form matrices AX, BX, CX, DX, V1X, and V2x bigax = bigpi; bigbx = bigw; bigcx = [ bigu(1,:) ; bigu(4,:) ; bigu(3,:) ]; bigdx = [ dyytr dyctr dyhtr ; ... dcytr dcctr dchtr ; ... dhytr dhctr dhhtr ]; dxeig = eig(bigdx); dxviol = 0; if max(abs(dxeig)) > 1 dxviol = 1; end bigv1x = sig^2; bigv2x1 = [ vyytr 0 0 ; ... vcytr vcctr 0 ; ... vhytr vhctr vhhtr ]; bigv2x = bigv2x1*bigv2x1’; % form matrices FX, GX, and QX bigfx = [ bigax zeros(2,3) ; zeros(3,2) bigdx ]; 21

biggx = [ bigcx eye(3) ]; bigqx = [ bigbx*bigv1x*bigbx’ zeros(2,3) ; zeros(3,2) bigv2x ]; % put data in deviation form trend = 1:bigt; ythat = log(yt) - log(eta)*trend’ - log(yss); cthat = log(ct) - log(eta)*trend’ - log(css); hthat = log(ht) - log(hss); dthat = [ ythat cthat hthat ]; % evaluate negative log likelihood xt = zeros(5,1); bigsig1 = inv(eye(25)-kron(bigfx,bigfx))*bigqx(:); bigsigt = reshape(bigsig1,5,5); llfn = (3*bigt/2)*log(2*pi); for t = 1:bigt ut = dthat(t,:)’ - biggx*xt; omegt = biggx*bigsigt*biggx’; omeginvt = inv(omegt); llfn = llfn + (1/2)*(log(det(omegt))+ut’*omeginvt*ut); bigkt = bigfx*bigsigt*biggx’*omeginvt; xt = bigfx*xt + bigkt*ut; bigsigt = bigqx + bigfx*bigsigt*bigfx’ ... - bigfx*bigsigt*biggx’*omeginvt*biggx*bigsigt*bigfx’; end % penalize constraint violations if dxviol llfn = llfn + 1e8; end if abs(imag(llfn)) > 0 llfn = real(llfn) + 1e8; end if abs(rhotr) > 1 llfn = llfn + 1e8; end

22

2.5. llfnse.m function llfn = llfnse(bigthet); % Uses the Kalman filter to evaluate the negative log likelihood % function for the real business cycle model with indivisible % labor. Parameters are untransformed to facilitate calculation % of standard errors. % % T H IS P R O G R A M WA S W R IT T E N FO R M AT L A B B Y P E T E R N . IR E L A N D % B O S T O N C O L L E G E , D E PA RT M E N T O F E C O N O M IC S 1 4 0 C O M M O N W E A LT H AV E N U E , % C H E S T N U T H IL L , M A 0 2 4 67 , irela n d p @ b c.e d u % % F IN A N C IA L S U P P O RT F R O M T H E N AT IO N A L S C IE N C E FO U N D AT IO N U N D E R G R A N T N O S . % S E S -9 9 8 57 6 3 A N D S E S -0 2 1 3 4 6 1 IS G R AT E F U L LY A C K N O W L E D G E D . % % C O P Y R IG H T (c) 20 0 3 B Y P E T E R N . IR E L A N D . R E % D IS T R IB U T IO N IS P E R M IT T E D FO R % E D U C AT IO N A L A N D R E S E A R C H P U R P O S E S , S O L O N G A S N O C H A N G E S A R E M A D E . A L L % C O P IE S M U S T B E P R O V ID E D F R E E O F C H A R G E A N D M U S T IN C L U D E T H IS C O P Y R IG H T N O T IC E .

% define variables and parameters global yt ct ht scalinv bigt = length(yt); bigthet = real(bigthet); bigthet = scalinv*bigthet; beta = 0.99; gamma = bigthet(1); theta = bigthet(2); eta = bigthet(3); delta = 0.025; a = bigthet(4); rho = bigthet(5); sig = bigthet(6); dyy = bigthet(7); dyc = bigthet(8); dyh = bigthet(9); dcy = bigthet(10); dcc = bigthet(11); 23

dch = bigthet(12); dhy = bigthet(13); dhc = bigthet(14); dhh = bigthet(15); vyy = bigthet(16); vcc = bigthet(17); vhh = bigthet(18); vyc = bigthet(19); vyh = bigthet(20); vch = bigthet(21); % compute steady-state values kappa = eta/beta - 1 + delta; lambda = eta - 1 + delta; hss = ((1-theta)/gamma)/(1-theta*lambda/kappa); yss = (a^(1/(1-theta)))*((theta/kappa)^(theta/(1-theta)))*hss; kss = (theta/kappa)*yss; iss = (theta*lambda/kappa)*yss; css = yss - iss; % compute K coefficients bigk11 = (eta-beta*(1-theta)*(1-delta))/(beta*eta*theta); bigk12 = (beta*eta*theta^2-eta+beta*(1-theta^2)*(1-delta)) ... /(beta*eta*theta^2); bigk22 = eta*theta/(eta-beta*(1-theta)*(1-delta)); % compute L coefficients bigl1 = (eta-beta*(1-delta))/(beta*eta*theta^2); bigl2 = rho*(eta-beta*(1-delta))/(eta-beta*(1-theta)*(1-delta)); % form S matrices bigs1 = (bigk22-bigk11)/bigk12; bigs2 = ((bigk22-bigk11)*bigl1-bigk12*bigl2)/(bigk12*(bigk11-rho)); bigs3 = bigk22; bigs4 = bigk12*bigl2/(bigk22-bigk11) ... + (bigk22-rho)*((bigk22-bigk11)*bigl1-bigk12*bigl2) ... / ((bigk22-bigk11)*(bigk11-rho)); bigs5 = [ 1 - ((1-theta)/theta)*bigs1 ; ... bigs1 + (kappa/(theta^2*lambda))*(theta-bigs1) ; ... 24

1 - (1/theta)*bigs1 ]; bigs6 = [ 1/theta - ((1-theta)/theta)*bigs2 ; ... bigs2 + (kappa/(theta^2*lambda))*(1-bigs2) ; ... 1/theta - (1/theta)*bigs2 ]; % form matrices PI, W, and U bigpi = [ bigs3 bigs4 ; 0 rho ]; bigw = [ 0 ; 1 ]; bigu = [ bigs5 bigs6 ; bigs1 bigs2 ]; % form matrices AX, BX, CX, DX, V1X, and V2x bigax = bigpi; bigbx = bigw; bigcx = [ bigu(1,:) ; bigu(4,:) ; bigu(3,:) ]; bigdx = [ dyy dyc dyh ; ... dcy dcc dch ; ... dhy dhc dhh ]; dxeig = eig(bigdx); dxviol = 0; if max(abs(dxeig)) > 1 dxviol = 1; end bigv1x = sig^2; bigv2x = [ vyy^2 vyc vyh ; ... vyc vcc^2 vch ; ... vyh vch vhh^2 ]; % form matrices FX, GX, and QX bigfx = [ bigax zeros(2,3) ; zeros(3,2) bigdx ]; biggx = [ bigcx eye(3) ]; bigqx = [ bigbx*bigv1x*bigbx’ zeros(2,3) ; zeros(3,2) bigv2x ]; % put data in deviation form trend = 1:bigt; ythat = log(yt) - log(eta)*trend’ - log(yss); cthat = log(ct) - log(eta)*trend’ - log(css); hthat = log(ht) - log(hss); dthat = [ ythat cthat hthat ]; % evaluate negative log likelihood 25

xt = zeros(5,1); bigsig1 = inv(eye(25)-kron(bigfx,bigfx))*bigqx(:); bigsigt = reshape(bigsig1,5,5); llfn = (3*bigt/2)*log(2*pi); for t = 1:bigt ut = dthat(t,:)’ - biggx*xt; omegt = biggx*bigsigt*biggx’; omeginvt = inv(omegt); llfn = llfn + (1/2)*(log(det(omegt))+ut’*omeginvt*ut); bigkt = bigfx*bigsigt*biggx’*omeginvt; xt = bigfx*xt + bigkt*ut; bigsigt = bigqx + bigfx*bigsigt*bigfx’ ... - bigfx*bigsigt*biggx’*omeginvt*biggx*bigsigt*bigfx’; end % penalize constraint violations if dxviol llfn = llfn + 1e8; end if abs(imag(llfn)) > 0 llfn = real(llfn) + 1e8; end if abs(rho) > 1 llfn = llfn + 1e8; end

26

3. Solution and Simulation 3.1. solv.m % solv.m: Solves the real business cycle model with indivisible labor. % % T H IS P R O G R A M WA S W R IT T E N FO R M AT L A B B Y P E T E R N . IR E L A N D % B O S T O N C O L L E G E , D E PA RT M E N T O F E C O N O M IC S 1 4 0 C O M M O N W E A LT H AV E N U E , % C H E S T N U T H IL L , M A 0 2 4 67 , irela n d p @ b c.e d u % % F IN A N C IA L S U P P O RT F R O M T H E N AT IO N A L S C IE N C E FO U N D AT IO N U N D E R G R A N T N O S . % S E S -9 9 8 57 6 3 A N D S E S -0 2 1 3 4 6 1 IS G R AT E F U L LY A C K N O W L E D G E D . % % C O P Y R IG H T (c) 20 0 3 B Y P E T E R N . IR E L A N D . R E % D IS T R IB U T IO N IS P E R M IT T E D FO R % E D U C AT IO N A L A N D R E S E A R C H P U R P O S E S , S O L O N G A S N O C H A N G E S A R E M A D E . A L L % C O P IE S M U S T B E P R O V ID E D F R E E O F C H A R G E A N D M U S T IN C L U D E T H IS C O P Y R IG H T N O T IC E .

% set parameter values beta = 0.99; gamma = 0.0045; theta = 0.20; eta = 1.0042; delta = 0.025; a = 6; rho = 0.95; sig = 0.007; % compute steady-state values kappa = eta/beta - 1 + delta; lambda = eta - 1 + delta; hss = ((1-theta)/gamma)/(1-theta*lambda/kappa); yss = (a^(1/(1-theta)))*((theta/kappa)^(theta/(1-theta)))*hss; kss = (theta/kappa)*yss; iss = (theta*lambda/kappa)*yss; css = yss - iss; % compute K coefficients bigk11 = (eta-beta*(1-theta)*(1-delta))/(beta*eta*theta); bigk12 = (beta*eta*theta^2-eta+beta*(1-theta^2)*(1-delta)) ... 27

/(beta*eta*theta^2); bigk22 = eta*theta/(eta-beta*(1-theta)*(1-delta)); % compute L coefficients bigl1 = (eta-beta*(1-delta))/(beta*eta*theta^2); bigl2 = rho*(eta-beta*(1-delta))/(eta-beta*(1-theta)*(1-delta)); % form S matrices bigs1 = (bigk22-bigk11)/bigk12; bigs2 = ((bigk22-bigk11)*bigl1-bigk12*bigl2)/(bigk12*(bigk11-rho)); bigs3 = bigk22; bigs4 = bigk12*bigl2/(bigk22-bigk11) ... + (bigk22-rho)*((bigk22-bigk11)*bigl1-bigk12*bigl2) ... / ((bigk22-bigk11)*(bigk11-rho)); bigs5 = [ 1 - ((1-theta)/theta)*bigs1 ; ... bigs1 + (kappa/(theta^2*lambda))*(theta-bigs1) ; ... 1 - (1/theta)*bigs1 ]; bigs6 = [ 1/theta - ((1-theta)/theta)*bigs2 ; ... bigs2 + (kappa/(theta^2*lambda))*(1-bigs2) ; ... 1/theta - (1/theta)*bigs2 ]; % form matrices PI, W, and U bigpi = [ bigs3 bigs4 ; 0 rho ]; bigw = [ 0 ; 1 ]; bigu = [ bigs5 bigs6 ; bigs1 bigs2 ];

28

3.2. check.m % check.m: Checks solutions found by solv.m for the real % business cycle model with indivisible labor % by making sure that the equilibrium conditions % hold in the steady state and after a technology % shock. If the solution is correct, the program % will return a 6x2 vector of zeros. % % T H IS P R O G R A M WA S W R IT T E N FO R M AT L A B B Y P E T E R N . IR E L A N D % B O S T O N C O L L E G E , D E PA RT M E N T O F E C O N O M IC S 1 4 0 C O M M O N W E A LT H AV E N U E , % C H E S T N U T H IL L , M A 0 2 4 67 , irela n d p @ b c.e d u % % F IN A N C IA L S U P P O RT F R O M T H E N AT IO N A L S C IE N C E FO U N D AT IO N U N D E R G R A N T N O S . % S E S -9 9 8 57 6 3 A N D S E S -0 2 1 3 4 6 1 IS G R AT E F U L LY A C K N O W L E D G E D . % % C O P Y R IG H T (c) 20 0 3 B Y P E T E R N . IR E L A N D . R E % D IS T R IB U T IO N IS P E R M IT T E D FO R % E D U C AT IO N A L A N D R E S E A R C H P U R P O S E S , S O L O N G A S N O C H A N G E S A R E M A D E . A L L % C O P IE S M U S T B E P R O V ID E D F R E E O F C H A R G E A N D M U S T IN C L U D E T H IS C O P Y R IG H T N O T IC E .

% check steady state conditions sscond = ones(6,1); sscond(1) = gamma*css*hss - (1-theta)*yss; sscond(2) = eta - beta*(theta*(yss/kss)+1-delta); sscond(3) = yss - css - iss; sscond(4) = yss - a*(kss^theta)*(hss^(1-theta)); sscond(5) = (eta-1+delta)*kss - iss; sscond(6) = a - a; % construct period t values after shock et = sig; st = bigw*et; ft = bigu*st; kt = st(1); at = st(2); yt = ft(1); it = ft(2); ht = ft(3); 29

ct = ft(4); % construct period t+1 values after shock stp = bigpi*st; ftp = bigu*stp; ktp = stp(1); atp = stp(2); ytp = ftp(1); itp = ftp(2); htp = ftp(3); ctp = ftp(4); % check equilibrium conditions eqcond = ones(6,1); eqcond(1) = ct + ht - yt; eqcond(2) = (eta/beta)*(ct-ctp) + (eta/beta-1+delta)*(ytp-ktp); eqcond(3) = (eta/beta-1+delta)*(yt-ct) + theta*(eta-1+delta)*(ct-it); eqcond(4) = yt - at - theta*kt - (1-theta)*ht; eqcond(5) = eta*ktp - (1-delta)*kt - (eta-1+delta)*it; eqcond(6) = atp - rho*at; % report results [ sscond eqcond ]

30

3.3. ksmooth.m % ksmooth.m: Uses the Kalman smoother, together with the un% transformed estimates found be estse.m, to produce a series % of smoothed estimates of the hybrid model’s shocks. % % The results are contained in the 5xT matrix xbigtvec, which % keeps track of the five state variables in the hybrid model’s % state equation. % % Also calculates corey, corec, coreh, the correlation % coefficients between the innovation to technology and the % innovations to each of the three residuals. % % T H IS P R O G R A M WA S W R IT T E N FO R M AT L A B B Y P E T E R N . IR E L A N D % B O S T O N C O L L E G E , D E PA RT M E N T O F E C O N O M IC S 1 4 0 C O M M O N W E A LT H AV E N U E , % C H E S T N U T H IL L , M A 0 2 4 67 , irela n d p @ b c.e d u % % F IN A N C IA L S U P P O RT F R O M T H E N AT IO N A L S C IE N C E FO U N D AT IO N U N D E R G R A N T N O S . % S E S -9 9 8 57 6 3 A N D S E S -0 2 1 3 4 6 1 IS G R AT E F U L LY A C K N O W L E D G E D . % % C O P Y R IG H T (c) 20 0 3 B Y P E T E R N . IR E L A N D . R E % D IS T R IB U T IO N IS P E R M IT T E D FO R % E D U C AT IO N A L A N D R E S E A R C H P U R P O S E S , S O L O N G A S N O C H A N G E S A R E M A D E . A L L % C O P IE S M U S T B E P R O V ID E D F R E E O F C H A R G E A N D M U S T IN C L U D E T H IS C O P Y R IG H T N O T IC E .

% compute steady-state values kappa = eta/beta - 1 + delta; lambda = eta - 1 + delta; hss = ((1-theta)/gamma)/(1-theta*lambda/kappa); yss = (a^(1/(1-theta)))*((theta/kappa)^(theta/(1-theta)))*hss; kss = (theta/kappa)*yss; iss = (theta*lambda/kappa)*yss; css = yss - iss; % compute K coefficients bigk11 = (eta-beta*(1-theta)*(1-delta))/(beta*eta*theta); bigk12 = (beta*eta*theta^2-eta+beta*(1-theta^2)*(1-delta)) ... /(beta*eta*theta^2); 31

bigk22 = eta*theta/(eta-beta*(1-theta)*(1-delta)); % compute L coefficients bigl1 = (eta-beta*(1-delta))/(beta*eta*theta^2); bigl2 = rho*(eta-beta*(1-delta))/(eta-beta*(1-theta)*(1-delta)); % form S matrices bigs1 = (bigk22-bigk11)/bigk12; bigs2 = ((bigk22-bigk11)*bigl1-bigk12*bigl2)/(bigk12*(bigk11-rho)); bigs3 = bigk22; bigs4 = bigk12*bigl2/(bigk22-bigk11) ... + (bigk22-rho)*((bigk22-bigk11)*bigl1-bigk12*bigl2) ... / ((bigk22-bigk11)*(bigk11-rho)); bigs5 = [ 1 - ((1-theta)/theta)*bigs1 ; ... bigs1 + (kappa/(theta^2*lambda))*(theta-bigs1) ; ... 1 - (1/theta)*bigs1 ]; bigs6 = [ 1/theta - ((1-theta)/theta)*bigs2 ; ... bigs2 + (kappa/(theta^2*lambda))*(1-bigs2) ; ... 1/theta - (1/theta)*bigs2 ]; % form matrices PI, W, and U bigpi = [ bigs3 bigs4 ; 0 rho ]; bigw = [ 0 ; 1 ]; bigu = [ bigs5 bigs6 ; bigs1 bigs2 ]; % form matrices AX, BX, CX, DX, V1X, and V2x bigax = bigpi; bigbx = bigw; bigcx = [ bigu(1,:) ; bigu(4,:) ; bigu(3,:) ]; bigdx = [ dyy dyc dyh ; ... dcy dcc dch ; ... dhy dhc dhh ]; dxeig = eig(bigdx); dxviol = 0; if max(abs(dxeig)) > 1 dxviol = 1; end bigv1x = sig^2; bigv2x = [ vyy^2 vyc vyh ; ... 32

vyc vcc^2 vch ; ... vyh vch vhh^2 ]; % form matrices FX, GX, and QX bigfx = [ bigax zeros(2,3) ; zeros(3,2) bigdx ]; biggx = [ bigcx eye(3) ]; bigqx = [ bigbx*bigv1x*bigbx’ zeros(2,3) ; zeros(3,2) bigv2x ]; % put data in deviation form bigt = length(yt); trend = 1:bigt; ythat = log(yt) - log(eta)*trend’ - log(yss); cthat = log(ct) - log(eta)*trend’ - log(css); hthat = log(ht) - log(hss); dthat = [ ythat cthat hthat ]; % run through the Kalman Filter xttvec = zeros(5,bigt); xtlvec = zeros(5,bigt); sigttvec = zeros(25,bigt); sigtlvec = zeros(25,bigt); xt = zeros(5,1); bigsigt = inv(eye(25)-kron(bigfx,bigfx))*bigqx(:); bigsigt = reshape(bigsigt,5,5); for t = 1:bigt xtlvec(:,t) = xt; sigtlvec(:,t) = bigsigt(:); omegt = biggx*bigsigt*biggx’; omeginvt = inv(omegt); ut = dthat(t,:)’ - biggx*xt; xtt = xt + bigsigt*biggx’*omeginvt*ut; bigsigtt = bigsigt - bigsigt*biggx’*omeginvt*biggx*bigsigt; xttvec(:,t) = xtt; sigttvec(:,t) = bigsigtt(:); xt = bigfx*xtt; bigsigt = bigqx + bigfx*bigsigtt*bigfx’; end % run through Kalman smoother 33

xbigtvec = zeros(5,bigt); xbigtvec(:,bigt) = xttvec(:,bigt); for j = 1:bigt-1 bigsigtt = sigttvec(:,bigt-j); bigsigtt = reshape(bigsigtt,5,5); bigsigtp = sigtlvec(:,bigt-j+1); bigsigtp = reshape(bigsigtp,5,5); bigjt = bigsigtt*bigfx’*inv(bigsigtp); xtt = xttvec(:,bigt-j); xtpbigt = xbigtvec(:,bigt-j+1); xtpt = xtlvec(:,bigt-j+1); xbigtvec(:,bigt-j)= xtt + bigjt*(xtpbigt-xtpt); end % calculate correlations between innovations etavec = xbigtvec(:,2:bigt)’ - xbigtvec(:,1:bigt-1)’*bigfx’; epsvec = etavec(:,2); xiyvec = etavec(:,3); xicvec = etavec(:,4); xihvec = etavec(:,5); corey = corrcoef(epsvec,xiyvec); corec = corrcoef(epsvec,xicvec); coreh = corrcoef(epsvec,xihvec);

34

4. Impulse Responses & Variance Decompostion 4.1. imp.m % imp.m: Computes impulse responses for the real business % cycle model with indivisible labor using the % solutions found by solv.m. % % T H IS P R O G R A M WA S W R IT T E N FO R M AT L A B B Y P E T E R N . IR E L A N D % B O S T O N C O L L E G E , D E PA RT M E N T O F E C O N O M IC S 1 4 0 C O M M O N W E A LT H AV E N U E , % C H E S T N U T H IL L , M A 0 2 4 67 , irela n d p @ b c.e d u % % F IN A N C IA L S U P P O RT F R O M T H E N AT IO N A L S C IE N C E FO U N D AT IO N U N D E R G R A N T N O S . % S E S -9 9 8 57 6 3 A N D S E S -0 2 1 3 4 6 1 IS G R AT E F U L LY A C K N O W L E D G E D . % % C O P Y R IG H T (c) 20 0 3 B Y P E T E R N . IR E L A N D . R E % D IS T R IB U T IO N IS P E R M IT T E D FO R % E D U C AT IO N A L A N D R E S E A R C H P U R P O S E S , S O L O N G A S N O C H A N G E S A R E M A D E . A L L % C O P IE S M U S T B E P R O V ID E D F R E E O F C H A R G E A N D M U S T IN C L U D E T H IS C O P Y R IG H T N O T IC E .

% technology shock st = zeros(50,2); ft = zeros(50,4); et = sig; st(2,:) = (bigw*et)’; ft(2,:) = (bigu*st(2,:)’)’; for t = 3:50; st(t,:) = (bigpi*st(t-1,:)’)’; ft(t,:) = (bigu*st(t,:)’)’; end % create output vectors st = 100*st; ft = 100*ft; kt = st(:,1); at = st(:,2); yt = ft(:,1); it = ft(:,2); ht = ft(:,3); 35

ct = ft(:,4);

36

4.2. vardec.m % vardec.m: Uses the parameter estimates and standard % errors found by estse.m to compute variance % decompositions and standard errors for the real % business cycle model with indivisible labor. % % Returns a 56x2 matrix vdec. The first column is: % % [ yvar cvar ivar hvar yvart cvart ivart hvart ]’ % % where yvar, cvar, ivar, and hvar are vectors of k-step % ahead forecast error variances in output, consumption, % investment, and hours worked, and k = 1,4,8,12,20,40, % and infinity. The variances are all expressed as % percentages. The vectors yvart, cvart, ivart, and % hvart are vectors of percentages due to the technology % shock. The second column gives the standard errors % for the corresponding elements of the first column. % % Formulas used to compute variance decompositions and % standard errors are contained in vardecfn.m. % % T H IS P R O G R A M WA S W R IT T E N FO R M AT L A B B Y P E T E R N . IR E L A N D % B O S T O N C O L L E G E , D E PA RT M E N T O F E C O N O M IC S 1 4 0 C O M M O N W E A LT H AV E N U E , % C H E S T N U T H IL L , M A 0 2 4 67 , irela n d p @ b c.e d u % % F IN A N C IA L S U P P O RT F R O M T H E N AT IO N A L S C IE N C E FO U N D AT IO N U N D E R G R A N T N O S . % S E S -9 9 8 57 6 3 A N D S E S -0 2 1 3 4 6 1 IS G R AT E F U L LY A C K N O W L E D G E D . % % C O P Y R IG H T (c) 20 0 3 B Y P E T E R N . IR E L A N D . R E % D IS T R IB U T IO N IS P E R M IT T E D FO R % E D U C AT IO N A L A N D R E S E A R C H P U R P O S E S , S O L O N G A S N O C H A N G E S A R E M A D E . A L L % C O P IE S M U S T B E P R O V ID E D F R E E O F C H A R G E A N D M U S T IN C L U D E T H IS C O P Y R IG H T N O T IC E .

% find covariance matrix of estimated parameters bighx = inv(hessmat); % compute variance decompositions 37

vdec1 = vardecfn(tstar); % compute standard errors epsmat = 1e-6*eye(21); gradmat = zeros(56,21); for i = 1:21 gradmat(:,i) = (vardecfn(tstar+epsmat(:,i))-vdec1)/1e-6; end vdec2 = gradmat*bighx*gradmat’; vdec2 = sqrt(diag(vdec2)); % collect results vdec = [ vdec1 vdec2 ];

38

4.3. vardecfn.m function vardecfn = vardecfn(bigthet); % For given parameter estimates, computes variance % decompositions for the real business cycle model % with indivisible labor. % % T H IS P R O G R A M WA S W R IT T E N FO R M AT L A B B Y P E T E R N . IR E L A N D % B O S T O N C O L L E G E , D E PA RT M E N T O F E C O N O M IC S 1 4 0 C O M M O N W E A LT H AV E N U E , % C H E S T N U T H IL L , M A 0 2 4 67 , irela n d p @ b c.e d u % % F IN A N C IA L S U P P O RT F R O M T H E N AT IO N A L S C IE N C E FO U N D AT IO N U N D E R G R A N T N O S . % S E S -9 9 8 57 6 3 A N D S E S -0 2 1 3 4 6 1 IS G R AT E F U L LY A C K N O W L E D G E D . % % C O P Y R IG H T (c) 20 0 3 B Y P E T E R N . IR E L A N D . R E % D IS T R IB U T IO N IS P E R M IT T E D FO R % E D U C AT IO N A L A N D R E S E A R C H P U R P O S E S , S O L O N G A S N O C H A N G E S A R E M A D E . A L L % C O P IE S M U S T B E P R O V ID E D F R E E O F C H A R G E A N D M U S T IN C L U D E T H IS C O P Y R IG H T N O T IC E .

% define parameters beta = 0.99; gamma = bigthet(1); theta = bigthet(2); eta = bigthet(3); delta = 0.025; a = bigthet(4); rho = bigthet(5); sig = bigthet(6); dyy = bigthet(7); dyc = bigthet(8); dyh = bigthet(9); dcy = bigthet(10); dcc = bigthet(11); dch = bigthet(12); dhy = bigthet(13); dhc = bigthet(14); dhh = bigthet(15); vyy = bigthet(16); 39

vcc = bigthet(17); vhh = bigthet(18); vyc = bigthet(19); vyh = bigthet(20); vch = bigthet(21); % compute steady-state values kappa = eta/beta - 1 + delta; lambda = eta - 1 + delta; hss = ((1-theta)/gamma)/(1-theta*lambda/kappa); yss = (a^(1/(1-theta)))*((theta/kappa)^(theta/(1-theta)))*hss; kss = (theta/kappa)*yss; iss = (theta*lambda/kappa)*yss; css = yss - iss; % compute K coefficients bigk11 = (eta-beta*(1-theta)*(1-delta))/(beta*eta*theta); bigk12 = (beta*eta*theta^2-eta+beta*(1-theta^2)*(1-delta)) ... /(beta*eta*theta^2); bigk22 = eta*theta/(eta-beta*(1-theta)*(1-delta)); % compute L coefficients bigl1 = (eta-beta*(1-delta))/(beta*eta*theta^2); bigl2 = rho*(eta-beta*(1-delta))/(eta-beta*(1-theta)*(1-delta)); % form S matrices bigs1 = (bigk22-bigk11)/bigk12; bigs2 = ((bigk22-bigk11)*bigl1-bigk12*bigl2)/(bigk12*(bigk11-rho)); bigs3 = bigk22; bigs4 = bigk12*bigl2/(bigk22-bigk11) ... + (bigk22-rho)*((bigk22-bigk11)*bigl1-bigk12*bigl2) ... / ((bigk22-bigk11)*(bigk11-rho)); bigs5 = [ 1 - ((1-theta)/theta)*bigs1 ; ... bigs1 + (kappa/(theta^2*lambda))*(theta-bigs1) ; ... 1 - (1/theta)*bigs1 ]; bigs6 = [ 1/theta - ((1-theta)/theta)*bigs2 ; ... bigs2 + (kappa/(theta^2*lambda))*(1-bigs2) ; ... 1/theta - (1/theta)*bigs2 ]; % form matrices PI, W, and U 40

bigpi = [ bigs3 bigs4 ; 0 rho ]; bigw = [ 0 ; 1 ]; bigu = [ bigs5 bigs6 ; bigs1 bigs2 ]; % form matrices AX, BX, CX, DX, V1X, and V2x bigax = bigpi; bigbx = bigw; bigcx = [ bigu(1,:) ; bigu(4,:) ; bigu(3,:) ]; bigdx = [ dyy dyc dyh ; ... dcy dcc dch ; ... dhy dhc dhh ]; dxeig = eig(bigdx); dxviol = 0; if max(abs(dxeig)) > 1 dxviol = 1; end bigv1x = sig^2; bigv2x = [ vyy^2 vyc vyh ; ... vyc vcc^2 vch ; ... vyh vch vhh^2 ]; % form matrices FX, GX, and QX bigfx = [ bigax zeros(2,3) ; zeros(3,2) bigdx ]; biggx = [ bigcx eye(3) ]; bigqx = [ bigbx*bigv1x*bigbx’ zeros(2,3) ; zeros(3,2) bigv2x ]; % calculate k-step ahead forecast error variances vardecfn = zeros(56,1); bigf2x = eye(5); bigg2x = zeros(4,5); bigg2x(1:2,:) = biggx(1:2,:); bigg2x(3,:) = (yss/iss)*biggx(1,:) - (css/iss)*biggx(2,:); bigg2x(4,:) = biggx(3,:); bigsigxk = zeros(5,5); for k = 1:40 bigsigxk = bigsigxk + bigf2x*bigqx*bigf2x’; bigsigyk = bigg2x*bigsigxk*bigg2x’; if k == 1 41

vardecfn(1) = bigsigyk(1,1); vardecfn(8) = bigsigyk(2,2); vardecfn(15) = bigsigyk(3,3); vardecfn(22) = bigsigyk(4,4); end if k == 4 vardecfn(2) = bigsigyk(1,1); vardecfn(9) = bigsigyk(2,2); vardecfn(16) = bigsigyk(3,3); vardecfn(23) = bigsigyk(4,4); end if k == 8 vardecfn(3) = bigsigyk(1,1); vardecfn(10) = bigsigyk(2,2); vardecfn(17) = bigsigyk(3,3); vardecfn(24) = bigsigyk(4,4); end if k == 12 vardecfn(4) = bigsigyk(1,1); vardecfn(11) = bigsigyk(2,2); vardecfn(18) = bigsigyk(3,3); vardecfn(25) = bigsigyk(4,4); end if k == 20 vardecfn(5) = bigsigyk(1,1); vardecfn(12) = bigsigyk(2,2); vardecfn(19) = bigsigyk(3,3); vardecfn(26) = bigsigyk(4,4); end if k == 40 vardecfn(6) = bigsigyk(1,1); vardecfn(13) = bigsigyk(2,2); vardecfn(20) = bigsigyk(3,3); vardecfn(27) = bigsigyk(4,4); end 42

bigf2x = bigf2x*bigfx; end bigsigx1 = inv(eye(25)-kron(bigfx,bigfx))*bigqx(:); bigsigx = reshape(bigsigx1,5,5); bigsigy = bigg2x*bigsigx*bigg2x’; vardecfn(7) = bigsigy(1,1); vardecfn(14) = bigsigy(2,2); vardecfn(21) = bigsigy(3,3); vardecfn(28) = bigsigy(4,4); % calculate fractions due to technology bigqxt = [ bigbx*bigv1x*bigbx’ zeros(2,3) ; zeros(3,5) ]; bigf2x = eye(5); bigsigxk = zeros(5,5); for k = 1:40 bigsigxk = bigsigxk + bigf2x*bigqxt*bigf2x’; bigsigyk = bigg2x*bigsigxk*bigg2x’; if k == 1 vardecfn(29) = bigsigyk(1,1)/vardecfn(1); vardecfn(36) = bigsigyk(2,2)/vardecfn(8); vardecfn(43) = bigsigyk(3,3)/vardecfn(15); vardecfn(50) = bigsigyk(4,4)/vardecfn(22); end if k == 4 vardecfn(30) = bigsigyk(1,1)/vardecfn(2); vardecfn(37) = bigsigyk(2,2)/vardecfn(9); vardecfn(44) = bigsigyk(3,3)/vardecfn(16); vardecfn(51) = bigsigyk(4,4)/vardecfn(23); end if k == 8 vardecfn(31) = bigsigyk(1,1)/vardecfn(3); vardecfn(38) = bigsigyk(2,2)/vardecfn(10); vardecfn(45) = bigsigyk(3,3)/vardecfn(17); vardecfn(52) = bigsigyk(4,4)/vardecfn(24); end if k == 12 43

vardecfn(32) = bigsigyk(1,1)/vardecfn(4); vardecfn(39) = bigsigyk(2,2)/vardecfn(11); vardecfn(46) = bigsigyk(3,3)/vardecfn(18); vardecfn(53) = bigsigyk(4,4)/vardecfn(25); end if k == 20 vardecfn(33) = bigsigyk(1,1)/vardecfn(5); vardecfn(40) = bigsigyk(2,2)/vardecfn(12); vardecfn(47) = bigsigyk(3,3)/vardecfn(19); vardecfn(54) = bigsigyk(4,4)/vardecfn(26); end if k == 40 vardecfn(34) = bigsigyk(1,1)/vardecfn(6); vardecfn(41) = bigsigyk(2,2)/vardecfn(13); vardecfn(48) = bigsigyk(3,3)/vardecfn(20); vardecfn(55) = bigsigyk(4,4)/vardecfn(27); end bigf2x = bigf2x*bigfx; end bigsigx1 = inv(eye(25)-kron(bigfx,bigfx))*bigqxt(:); bigsigx = reshape(bigsigx1,5,5); bigsigy = bigg2x*bigsigx*bigg2x’; vardecfn(35) = bigsigy(1,1)/vardecfn(7); vardecfn(42) = bigsigy(2,2)/vardecfn(14); vardecfn(49) = bigsigy(3,3)/vardecfn(21); vardecfn(56) = bigsigy(4,4)/vardecfn(28); % convert to percentages vardecfn = 100*vardecfn;

44

5. Stability Tests 5.1. stabtest.m % stabtest.m: Performs parameter stability tests for % the real business cycle model with indivisible labor. % Breakpoint occurs at 1973:1. % % Returns: % % tstar1 = estimated parameters, 1948:1-1972:4 % sevec1 = standard errors, 1948:1-1972:4 % tstar2 = estimated parameters, 1973:1-2002:2 % sevec2 = standard errors, 1973:1-2002:2 % lrstat = LR statistic for stability of all parameters % wstat = Wald statistic for stability of all parameters % wstatq = Wald statistic for stability of structural % parameters % wstatx = Wald statistic for stability of structural % parameters, excluding average productivity growth % wstatn = Wald statistic for stability of other % parameters % % T H IS P R O G R A M WA S W R IT T E N FO R M AT L A B B Y P E T E R N . IR E L A N D % B O S T O N C O L L E G E , D E PA RT M E N T O F E C O N O M IC S 1 4 0 C O M M O N W E A LT H AV E N U E , % C H E S T N U T H IL L , M A 0 2 4 67 , irela n d p @ b c.e d u % % F IN A N C IA L S U P P O RT F R O M T H E N AT IO N A L S C IE N C E FO U N D AT IO N U N D E R G R A N T N O S . % S E S -9 9 8 57 6 3 A N D S E S -0 2 1 3 4 6 1 IS G R AT E F U L LY A C K N O W L E D G E D . % % C O P Y R IG H T (c) 20 0 3 B Y P E T E R N . IR E L A N D . R E % D IS T R IB U T IO N IS P E R M IT T E D FO R % E D U C AT IO N A L A N D R E S E A R C H P U R P O S E S , S O L O N G A S N O C H A N G E S A R E M A D E . A L L % C O P IE S M U S T B E P R O V ID E D F R E E O F C H A R G E A N D M U S T IN C L U D E T H IS C O P Y R IG H T N O T IC E .

% load data global yt ct ht scalinv load ych.dat; 45

yt = ych(:,1); ct = ych(:,2); ht = ych(:,3); % define pre-1973 subsample yt = yt(1:100); ct = ct(1:100); ht = ht(1:100); % set starting values for pre-1973 subsample bettr = sqrt(0.99/(1-0.99)); gamtr = 0.0045; thettr = sqrt(0.20/(1-0.20)); etatr = 0.0051; delttr = sqrt(0.025/(1-0.025)); atr = 6; rhotr = 0.9964; sigtr = 0.0054; dyytr = 1.1883; dyctr = 0.3612; dyhtr = -0.4075; dcytr = 0.0875; dcctr = 1.0024; dchtr = -0.0738; dhytr = 0.6527; dhctr = 0.3161; dhhtr = 0.3644; vyytr = 0.0090; vcytr = 0.0061; vcctr = -0.0064; vhytr = 0.0029; vhctr = -0.0011; vhhtr = -0.0031; % maximize likelihood for pre-1973 subsample bigtheto = [ gamtr thettr etatr atr rhotr sigtr ... dyytr dyctr dyhtr ... dcytr dcctr dchtr ... 46

dhytr dhctr dhhtr ... vyytr ... vcytr vcctr ... vhytr vhctr vhhtr ]’; options(1) = 1; options(14) = 10000; thetstar = fminu(’llfn’,bigtheto,options); % find standard errors for pre-1973 subsample thetstar = real(thetstar); beta = bettr^2/(1+bettr^2); gamma = abs(thetstar(1)); theta = thetstar(2)^2/(1+thetstar(2)^2); eta = 1 + abs(thetstar(3)); delta = delttr^2/(1+delttr^2); a = abs(thetstar(4)); rho = thetstar(5); sig = abs(thetstar(6)); dyy = thetstar(7); dyc = thetstar(8); dyh = thetstar(9); dcy = thetstar(10); dcc = thetstar(11); dch = thetstar(12); dhy = thetstar(13); dhc = thetstar(14); dhh = thetstar(15); cholv = [ thetstar(16) 0 0 ; ... thetstar(17:18)’ 0 ; ... thetstar(19:21)’ ]; bigv = cholv*cholv’; vyy = sqrt(bigv(1,1)); vcc = sqrt(bigv(2,2)); vhh = sqrt(bigv(3,3)); vyc = bigv(1,2); vyh = bigv(1,3); 47

vch = bigv(2,3); tstar1 = [ gamma theta eta a rho sig ... dyy dyc dyh dcy dcc dch dhy dhc dhh ... vyy vcc vhh vyc vyh vch ]’; scalvec = [ ones(19,1) ; 10 ; 10 ]; scalinv = inv(diag(scalvec)); tstar1s = diag(scalvec)*tstar1; fstar1 = llfnse(tstar1s); eee = 1e-6; epsmat = eee*eye(21); hessvec = zeros(21,1); for i = 1:21 hessvec(i) = llfnse(tstar1s+epsmat(:,i)); end hessmat = zeros(21,21); for i = 1:21 for j = 1:21 hessmat(i,j) = (llfnse(tstar1s+epsmat(:,i)+epsmat(:,j)) ... -hessvec(i)-hessvec(j)+fstar1)/eee^2; end end bighx1 = scalinv*inv(hessmat)*scalinv’; sevec1 = sqrt(diag(bighx1)); % reload data load ych.dat; yt = ych(:,1); ct = ych(:,2); ht = ych(:,3); % define post-1973 subsample yt = yt(101:218); ct = ct(101:218); ht = ht(101:218); % set starting values for post-1973 subsample bettr = sqrt(0.99/(1-0.99)); gamtr = 0.0045; 48

thettr = sqrt(0.20/(1-0.20)); etatr = 0.0051; delttr = sqrt(0.025/(1-0.025)); atr = 6; rhotr = 0.9991; sigtr = 0.0048; dyytr = 1.2245; dyctr = 0.7933; dyhtr = -0.4947; dcytr = 0.0285; dcctr = 1.0665; dchtr = -0.0560; dhytr = 0.3008; dhctr = 0.9692; dhhtr = 0.3935; vyytr = 0.0057; vcytr = 0.0034; vcctr = -0.0041; vhytr = 0.0012; vhctr = -0.0015; vhhtr = 0.0000; % maximize likelihood for post-1973 subsample bigtheto = [ gamtr thettr etatr atr rhotr sigtr ... dyytr dyctr dyhtr ... dcytr dcctr dchtr ... dhytr dhctr dhhtr ... vyytr ... vcytr vcctr ... vhytr vhctr vhhtr ]’; options(1) = 1; options(14) = 10000; thetstar = fminu(’llfn2’,bigtheto,options); % find standard errors for post-1973 subsample thetstar = real(thetstar); beta = bettr^2/(1+bettr^2); 49

gamma = abs(thetstar(1)); theta = thetstar(2)^2/(1+thetstar(2)^2); eta = 1 + abs(thetstar(3)); delta = delttr^2/(1+delttr^2); a = abs(thetstar(4)); rho = thetstar(5); sig = abs(thetstar(6)); dyy = thetstar(7); dyc = thetstar(8); dyh = thetstar(9); dcy = thetstar(10); dcc = thetstar(11); dch = thetstar(12); dhy = thetstar(13); dhc = thetstar(14); dhh = thetstar(15); cholv = [ thetstar(16) 0 0 ; ... thetstar(17:18)’ 0 ; ... thetstar(19:21)’ ]; bigv = cholv*cholv’; vyy = sqrt(bigv(1,1)); vcc = sqrt(bigv(2,2)); vhh = sqrt(bigv(3,3)); vyc = bigv(1,2); vyh = bigv(1,3); vch = bigv(2,3); tstar2 = [ gamma theta eta a rho sig ... dyy dyc dyh dcy dcc dch dhy dhc dhh ... vyy vcc vhh vyc vyh vch ]’; scalvec = ones(21,1); scalinv = inv(diag(scalvec)); tstar2s = diag(scalvec)*tstar2; fstar2 = llfnse2(tstar2); eee = 1e-6; epsmat = eee*eye(21); 50

hessvec = zeros(21,1); for i = 1:21 hessvec(i) = llfnse2(tstar2s+epsmat(:,i)); end hessmat = zeros(21,21); for i = 1:21 for j = 1:21 hessmat(i,j) = (llfnse2(tstar2s+epsmat(:,i)+epsmat(:,j)) ... -hessvec(i)-hessvec(j)+fstar2)/eee^2; end end bighx2 = scalinv*inv(hessmat)*scalinv’; sevec2 = sqrt(diag(bighx2)); % form test statistics lrstat = -2*(fstar1+fstar2+2323.5501); wstat = (tstar1-tstar2)’*inv(bighx1+bighx2)*(tstar1-tstar2); bighxq1 = bighx1(1:6,1:6); bighxq2 = bighx2(1:6,1:6); wstatq = (tstar1(1:6)-tstar2(1:6))’ ... *inv(bighxq1+bighxq2)*(tstar1(1:6)-tstar2(1:6)); bighxx1 = [ bighx1(1:2,1:2) bighx1(1:2,4:6) ; ... bighx1(4:6,1:2) bighx1(4:6,4:6) ]; bighxx2 = [ bighx2(1:2,1:2) bighx2(1:2,4:6) ; ... bighx2(4:6,1:2) bighx2(4:6,4:6) ]; tstarx1 = [ tstar1(1:2) ; tstar1(4:6) ]; tstarx2 = [ tstar2(1:2) ; tstar2(4:6) ]; wstatx = (tstarx1-tstarx2)’ ... *inv(bighxx1+bighxx2)*(tstarx1-tstarx2); bighxn1 = bighx1(7:21,7:21); bighxn2 = bighx2(7:21,7:21); wstatn = (tstar1(7:21)-tstar2(7:21))’ ... *inv(bighxn1+bighxn2)*(tstar1(7:21)-tstar2(7:21));

51

6. Forecasting 6.1. fork.m % fork.m: Generates k-step-ahead forecasts from four competing models: % % 1) The real business cycle model with indivisible labor, % where the matrices D and V are unconstrained. % % 2) The real business cycle model with indivisible labor, % where the matrices D and V are diagonal. % % 3) An unconstrained VAR(1). % % 4) An unconstrained VAR(2). % % Forecasts start using the model estimated over the sample % period 1948:1-1984:4 and end using the model estimated % over the sample period 1948:1-2002:2. % % Requires tstarmat and tstrmatd, the matrices of estimated % coefficients for the two versions of the business cycle % model generated by estseq.m and estseqd,m. Uses the % function forkfn.m to generate the forecasts from the % business cycle models. % % Returns: % % forkmatm = forecasts of [log(Y) log(C) log(I) log(H)] % from the unconstrained business cycle model. % % forkmatd = forecasts of [log(Y) log(C) log(I) log(H)] % from the constrained business cycle model. % % forkmatu = forecasts of [log(Y) log(C) log(I) log(H)] % from the VAR(1). 52

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %

forkmat2 = forecasts of [log(Y) log(C) log(I) log(H)] from the VAR(2). rmsem = root mean squared errors, unconstrained business cycle model rmsed = root mean squared errors, constrained business cycle model rmseu = root mean squared errors, VAR(1) rmse2 = root mean squared errors, VAR(2) ldm = mean loss differentials, unconstrained business cycle model, using the VAR(1) as a benchmark ldsem = standard error of ldm ldd = mean loss differentials, constrained business cycle model, using the VAR(1) as a benchmark ldsed = standard error of ldd ldm2 = mean loss differentials, unconstrained business cycle model, using the VAR(2) as a benchmark ldsem2 = standard error of ldm2 ldd2 = mean loss differentials, constrained business cycle model, using the VAR(2) as a benchmark ldsed2 = standard error of ldd2

% T H IS P R O G R A M WA S W R IT T E N FO R M AT L A B B Y P E T E R N . IR E L A N D

53

% B O S T O N C O L L E G E , D E PA RT M E N T O F E C O N O M IC S 1 4 0 C O M M O N W E A LT H AV E N U E , % C H E S T N U T H IL L , M A 0 2 4 67 , irela n d p @ b c.e d u % % F IN A N C IA L S U P P O RT F R O M T H E N AT IO N A L S C IE N C E FO U N D AT IO N U N D E R G R A N T N O S . % S E S -9 9 8 57 6 3 A N D S E S -0 2 1 3 4 6 1 IS G R AT E F U L LY A C K N O W L E D G E D . % % C O P Y R IG H T (c) 20 0 3 B Y P E T E R N . IR E L A N D . R E % D IS T R IB U T IO N IS P E R M IT T E D FO R % E D U C AT IO N A L A N D R E S E A R C H P U R P O S E S , S O L O N G A S N O C H A N G E S A R E M A D E . A L L % C O P IE S M U S T B E P R O V ID E D F R E E O F C H A R G E A N D M U S T IN C L U D E T H IS C O P Y R IG H T N O T IC E .

% load data and parameter estimates global yt ct ht load ych.dat; load tstarmat; load tstrmatd; % set forecast horizon and set up output matrices k = 4; forkmatm = zeros(71,4); forkmatd = zeros(71,4); forkmatu = zeros(71,4); forkmat2 = zeros(71,4); % generate forecasts from unconstrained business cycle model for j = 1:71 yt = ych(1:147+j,1); ct = ych(1:147+j,2); ht = ych(1:147+j,3); tstar = [ tstarmat(:,j) ; k ]; forkmatm(j,:) = forkfn(tstar); end % generate forecasts from constrained business cycle model for j = 1:71 yt = ych(1:147+j,1); ct = ych(1:147+j,2); ht = ych(1:147+j,3); tstar = [ tstrmatd(:,j) ; k ]; forkmatd(j,:) = forkfn(tstar); 54

end % generate forecasts from VAR(1) for j = 1:71 lyt = log(ych(1:147+j,1)); lct = log(ych(1:147+j,2)); lht = log(ych(1:147+j,3)); trend = 1:147+j; xx = [ ones(147+j,1) trend’ ]; bety = inv(xx’*xx)*xx’*lyt; betc = inv(xx’*xx)*xx’*lct; beth = inv(xx’*xx)*xx’*lht; ybar = bety(1); ytrend = bety(2); cbar = betc(1); ctrend = betc(2); hbar = beth(1); htrend = beth(2); dlyt = lyt - ybar - ytrend*trend’; dlct = lct - cbar - ctrend*trend’; dlht = lht - hbar - htrend*trend’; xxx = [ dlyt(1:146+j) dlct(1:146+j) dlht(1:146+j) ]; yyy = [ dlyt(2:147+j) dlct(2:147+j) dlht(2:147+j) ]; bigrho = (inv(xxx’*xxx)*xxx’*yyy)’; bigppp = bigrho; dlychtk = (bigrho^k)*yyy(146+j,:)’; dlytk = dlychtk(1); dlctk = dlychtk(2); dlhtk = dlychtk(3); lytk = dlytk + ybar + ytrend*(147+j+k); lctk = dlctk + cbar + ctrend*(147+j+k); litk = log(exp(lytk)-exp(lctk)); lhtk = dlhtk + hbar + htrend*(147+j+k); forkmatu(j,:) = [ lytk lctk litk lhtk ]; end % generate forecasts from VAR(2) 55

for j = 1:71 lyt = log(ych(1:147+j,1)); lct = log(ych(1:147+j,2)); lht = log(ych(1:147+j,3)); trend = 1:147+j; xx = [ ones(147+j,1) trend’ ]; bety = inv(xx’*xx)*xx’*lyt; betc = inv(xx’*xx)*xx’*lct; beth = inv(xx’*xx)*xx’*lht; ybar = bety(1); ytrend = bety(2); cbar = betc(1); ctrend = betc(2); hbar = beth(1); htrend = beth(2); dlyt = lyt - ybar - ytrend*trend’; dlct = lct - cbar - ctrend*trend’; dlht = lht - hbar - htrend*trend’; xxx = [ dlyt(2:146+j) dlct(2:146+j) dlht(2:146+j) ... dlyt(1:145+j) dlct(1:145+j) dlht(1:145+j) ]; yyy = [ dlyt(3:147+j) dlct(3:147+j) dlht(3:147+j) ... dlyt(2:146+j) dlct(2:146+j) dlht(2:146+j) ]; bigrho = (inv(xxx’*xxx)*xxx’*yyy)’; dlychtk = (bigrho^k)*yyy(145+j,:)’; dlytk = dlychtk(1); dlctk = dlychtk(2); dlhtk = dlychtk(3); lytk = dlytk + ybar + ytrend*(147+j+k); lctk = dlctk + cbar + ctrend*(147+j+k); litk = log(exp(lytk)-exp(lctk)); lhtk = dlhtk + hbar + htrend*(147+j+k); forkmat2(j,:) = [ lytk lctk litk lhtk ]; end % compute forecast errors lyt = log(ych(148+k:218,1)); 56

lct = log(ych(148+k:218,2)); lit = log(exp(lyt)-exp(lct)); lht = log(ych(148+k:218,3)); ldt = [ lyt lct lit lht ]; forerrm = forkmatm(1:71-k,:) - ldt; forerrd = forkmatd(1:71-k,:) - ldt; forerru = forkmatu(1:71-k,:) - ldt; forerr2 = forkmat2(1:71-k,:) - ldt; ldiffm = forerru.^2 - forerrm.^2; ldiffmy = ldiffm(:,1); ldiffmc = ldiffm(:,2); ldiffmi = ldiffm(:,3); ldiffmh = ldiffm(:,4); ldiffd = forerru.^2 - forerrd.^2; ldiffdy = ldiffd(:,1); ldiffdc = ldiffd(:,2); ldiffdi = ldiffd(:,3); ldiffdh = ldiffd(:,4); ldiffm2 = forerr2.^2 - forerrm.^2; ldiffm2y = ldiffm2(:,1); ldiffm2c = ldiffm2(:,2); ldiffm2i = ldiffm2(:,3); ldiffm2h = ldiffm2(:,4); ldiffd2 = forerr2.^2 - forerrd.^2; ldiffd2y = ldiffd2(:,1); ldiffd2c = ldiffd2(:,2); ldiffd2i = ldiffd2(:,3); ldiffd2h = ldiffd2(:,4); % compute forecast performance statistics rmsem = 100*sqrt((1/(71-k))*diag(forerrm’*forerrm)); rmsed = 100*sqrt((1/(71-k))*diag(forerrd’*forerrd)); rmseu = 100*sqrt((1/(71-k))*diag(forerru’*forerru)); rmse2 = 100*sqrt((1/(71-k))*diag(forerr2’*forerr2)); ldm = mean(ldiffm)’; ldvarmy = (1/(71-k))*ldiffmy’*ldiffmy; 57

ldvarmc = (1/(71-k))*ldiffmc’*ldiffmc; ldvarmi = (1/(71-k))*ldiffmi’*ldiffmi; ldvarmh = (1/(71-k))*ldiffmh’*ldiffmh; for l = 1:k-1 ldvarmy = ldvarmy + (2/(71-k))*ldiffmy(l+1:71-k)’*ldiffmy(1:71-k-l); ldvarmc = ldvarmc + (2/(71-k))*ldiffmc(l+1:71-k)’*ldiffmc(1:71-k-l); ldvarmi = ldvarmi + (2/(71-k))*ldiffmi(l+1:71-k)’*ldiffmi(1:71-k-l); ldvarmh = ldvarmh + (2/(71-k))*ldiffmh(l+1:71-k)’*ldiffmh(1:71-k-l); end ldvarm = (1/(71-k))*[ ldvarmy ; ldvarmc ; ldvarmi ; ldvarmh ]; ldsem = sqrt(ldvarm); ldd = mean(ldiffd)’; ldvardy = (1/(71-k))*ldiffdy’*ldiffdy; ldvardc = (1/(71-k))*ldiffdc’*ldiffdc; ldvardi = (1/(71-k))*ldiffdi’*ldiffdi; ldvardh = (1/(71-k))*ldiffdh’*ldiffdh; for l = 1:k-1 ldvardy = ldvardy + (2/(71-k))*ldiffdy(l+1:71-k)’*ldiffdy(1:71-k-l); ldvardc = ldvardc + (2/(71-k))*ldiffdc(l+1:71-k)’*ldiffdc(1:71-k-l); ldvardi = ldvardi + (2/(71-k))*ldiffdi(l+1:71-k)’*ldiffdi(1:71-k-l); ldvardh = ldvardh + (2/(71-k))*ldiffdh(l+1:71-k)’*ldiffdh(1:71-k-l); end ldvard = (1/(71-k))*[ ldvardy ; ldvardc ; ldvardi ; ldvardh ]; ldsed = sqrt(ldvard); ldm2 = mean(ldiffm2)’; ldvarm2y = (1/(71-k))*ldiffm2y’*ldiffm2y; ldvarm2c = (1/(71-k))*ldiffm2c’*ldiffm2c; ldvarm2i = (1/(71-k))*ldiffm2i’*ldiffm2i; ldvarm2h = (1/(71-k))*ldiffm2h’*ldiffm2h; for l = 1:k-1 ldvarm2y = ldvarm2y + (2/(71-k))*ldiffm2y(l+1:71-k)’*ldiffm2y(1:71-k-l); ldvarm2c = ldvarm2c + (2/(71-k))*ldiffm2c(l+1:71-k)’*ldiffm2c(1:71-k-l); ldvarm2i = ldvarm2i + (2/(71-k))*ldiffm2i(l+1:71-k)’*ldiffm2i(1:71-k-l); ldvarm2h = ldvarm2h + (2/(71-k))*ldiffm2h(l+1:71-k)’*ldiffm2h(1:71-k-l); end 58

ldvarm2 = (1/(71-k))*[ ldvarm2y ; ldvarm2c ; ldvarm2i ; ldvarm2h ]; ldsem2 = sqrt(ldvarm2); ldd2 = mean(ldiffd2)’; ldvard2y = (1/(71-k))*ldiffd2y’*ldiffd2y; ldvard2c = (1/(71-k))*ldiffd2c’*ldiffd2c; ldvard2i = (1/(71-k))*ldiffd2i’*ldiffd2i; ldvard2h = (1/(71-k))*ldiffd2h’*ldiffd2h; for l = 1:k-1 ldvard2y = ldvard2y + (2/(71-k))*ldiffd2y(l+1:71-k)’*ldiffd2y(1:71-k-l); ldvard2c = ldvard2c + (2/(71-k))*ldiffd2c(l+1:71-k)’*ldiffd2c(1:71-k-l); ldvard2i = ldvard2i + (2/(71-k))*ldiffd2i(l+1:71-k)’*ldiffd2i(1:71-k-l); ldvard2h = ldvard2h + (2/(71-k))*ldiffd2h(l+1:71-k)’*ldiffd2h(1:71-k-l); end ldvard2 = (1/(71-k))*[ ldvard2y ; ldvard2c ; ldvard2i ; ldvard2h ]; ldsed2 = sqrt(ldvard2);

59

6.2. forkfn.m function fork = forkfn(bigthet); % Uses the Kalman filter to generate k-step-ahead forecasts for output, % consumption, investment, and hours worked from the real business % cycle model with indivisible labor % % T H IS P R O G R A M WA S W R IT T E N FO R M AT L A B B Y P E T E R N . IR E L A N D % B O S T O N C O L L E G E , D E PA RT M E N T O F E C O N O M IC S 1 4 0 C O M M O N W E A LT H AV E N U E , % C H E S T N U T H IL L , M A 0 2 4 67 , irela n d p @ b c.e d u % % F IN A N C IA L S U P P O RT F R O M T H E N AT IO N A L S C IE N C E FO U N D AT IO N U N D E R G R A N T N O S . % S E S -9 9 8 57 6 3 A N D S E S -0 2 1 3 4 6 1 IS G R AT E F U L LY A C K N O W L E D G E D . % % C O P Y R IG H T (c) 20 0 3 B Y P E T E R N . IR E L A N D . R E % D IS T R IB U T IO N IS P E R M IT T E D FO R % E D U C AT IO N A L A N D R E S E A R C H P U R P O S E S , S O L O N G A S N O C H A N G E S A R E M A D E . A L L % C O P IE S M U S T B E P R O V ID E D F R E E O F C H A R G E A N D M U S T IN C L U D E T H IS C O P Y R IG H T N O T IC E .

% define variables and parameters global yt ct ht bigt = length(yt); bigthet = real(bigthet); beta = 0.99; gamma = bigthet(1); theta = bigthet(2); eta = bigthet(3); delta = 0.025; a = bigthet(4); rho = bigthet(5); sig = bigthet(6); dyy = bigthet(7); dyc = bigthet(8); dyh = bigthet(9); dcy = bigthet(10); dcc = bigthet(11); dch = bigthet(12); dhy = bigthet(13); 60

dhc = bigthet(14); dhh = bigthet(15); vyy = bigthet(16); vcc = bigthet(17); vhh = bigthet(18); vyc = bigthet(19); vyh = bigthet(20); vch = bigthet(21); k = bigthet(22); % compute steady-state values kappa = eta/beta - 1 + delta; lambda = eta - 1 + delta; hss = ((1-theta)/gamma)/(1-theta*lambda/kappa); yss = (a^(1/(1-theta)))*((theta/kappa)^(theta/(1-theta)))*hss; kss = (theta/kappa)*yss; iss = (theta*lambda/kappa)*yss; css = yss - iss; % compute K coefficients bigk11 = (eta-beta*(1-theta)*(1-delta))/(beta*eta*theta); bigk12 = (beta*eta*theta^2-eta+beta*(1-theta^2)*(1-delta)) ... /(beta*eta*theta^2); bigk22 = eta*theta/(eta-beta*(1-theta)*(1-delta)); % compute L coefficients bigl1 = (eta-beta*(1-delta))/(beta*eta*theta^2); bigl2 = rho*(eta-beta*(1-delta))/(eta-beta*(1-theta)*(1-delta)); % form S matrices bigs1 = (bigk22-bigk11)/bigk12; bigs2 = ((bigk22-bigk11)*bigl1-bigk12*bigl2)/(bigk12*(bigk11-rho)); bigs3 = bigk22; bigs4 = bigk12*bigl2/(bigk22-bigk11) ... + (bigk22-rho)*((bigk22-bigk11)*bigl1-bigk12*bigl2) ... / ((bigk22-bigk11)*(bigk11-rho)); bigs5 = [ 1 - ((1-theta)/theta)*bigs1 ; ... bigs1 + (kappa/(theta^2*lambda))*(theta-bigs1) ; ... 1 - (1/theta)*bigs1 ]; 61

bigs6 = [ 1/theta - ((1-theta)/theta)*bigs2 ; ... bigs2 + (kappa/(theta^2*lambda))*(1-bigs2) ; ... 1/theta - (1/theta)*bigs2 ]; % form matrices PI, W, and U bigpi = [ bigs3 bigs4 ; 0 rho ]; bigw = [ 0 ; 1 ]; bigu = [ bigs5 bigs6 ; bigs1 bigs2 ]; % form matrices AX, BX, CX, DX, V1X, and V2x bigax = bigpi; bigbx = bigw; bigcx = [ bigu(1,:) ; bigu(4,:) ; bigu(3,:) ]; bigdx = [ dyy dyc dyh ; ... dcy dcc dch ; ... dhy dhc dhh ]; dxeig = eig(bigdx); dxviol = 0; if max(abs(dxeig)) > 1 dxviol = 1; end bigv1x = sig^2; bigv2x = [ vyy^2 vyc vyh ; ... vyc vcc^2 vch ; ... vyh vch vhh^2 ]; % form matrices FX, GX, and QX bigfx = [ bigax zeros(2,3) ; zeros(3,2) bigdx ]; biggx = [ bigcx eye(3) ]; bigqx = [ bigbx*bigv1x*bigbx’ zeros(2,3) ; zeros(3,2) bigv2x ]; % put data in deviation form trend = 1:bigt; ythat = log(yt) - log(eta)*trend’ - log(yss); cthat = log(ct) - log(eta)*trend’ - log(css); hthat = log(ht) - log(hss); dthat = [ ythat cthat hthat ]; % generate forecasts of transformed variables xt = zeros(5,1); 62

bigsig1 = inv(eye(25)-kron(bigfx,bigfx))*bigqx(:); bigsigt = reshape(bigsig1,5,5); for t = 1:bigt ut = dthat(t,:)’ - biggx*xt; omegt = biggx*bigsigt*biggx’; omeginvt = inv(omegt); bigkt = bigfx*bigsigt*biggx’*omeginvt; xt = bigfx*xt + bigkt*ut; bigsigt = bigqx + bigfx*bigsigt*bigfx’ ... - bigfx*bigsigt*biggx’*omeginvt*biggx*bigsigt*bigfx’; end ytk = biggx*(bigfx^(k-1))*xt; % generate forecasts of untransformed variables bigytk = ytk(1) + (bigt+k)*log(eta) + log(yss); bigctk = ytk(2) + (bigt+k)*log(eta) + log(css); bigitk = log(exp(bigytk)-exp(bigctk)); bightk = ytk(3) + log(hss); fork = [ bigytk bigctk bigitk bightk ];

63