%------------------------------------------------------------------------------- % Article: "Persistence, Excesss Volatility, and Volatility Clusters in Inflation % Federal Reserve Bank of St Louis Review % Author: Michael Owyang % Issue: November/December 2001 % This file contains the data and calculations for: % Tables 1-3 %------------------------------------------------------------------------------- %------------------------------------------------------------------------------- % % Function that generates Zt (10/8/99) % Zt is the preference state indicator % Also generates the conditional probability of each state % Does this generate the transition probability? % Keep in mind that conditional on beta, y, and the hyperparameters, % the function that generates Zt is independent of the transition % equation and its components. % %------------------------------------------------------------------------------- function [ZTT,ZTT2]=genzt3d(param,trparamq,trparamp,betatt,STT); load data48; yy=dat48(206:600,:); T=size(yy,1); % Extract the parameters from the input vector sigett=param(1,1); sigvtt=param(2,1); eta2tt=param(3,1); eta1tt=param(4,1); alfa3tt=param(5,1); alfa2tt=param(6,1); alfa1tt=param(7,1); ktt=param(8,1); q33tt=trparamq(3,3); q22tt=trparamq(2,2); q11tt=trparamq(1,1); q31tt=trparamq(3,1); q32tt=trparamq(3,2); q13tt=trparamq(1,3); q12tt=trparamq(1,2); q23tt=trparamq(2,3); q21tt=trparamq(2,1); p22tt=trparamp(2,2); p11tt=trparamp(1,1); p12tt=trparamp(1,2); p21tt=trparamp(2,1); trmat=[q11tt,q12tt ,q13tt ,q21tt ,q22tt, q23tt , q31tt , q32tt ,q33tt]'; % transition matrix % Initial condititons q33ss=(2-q22tt-q11tt)/(3-q11tt-q22tt-q33tt); q22ss=(2-q11tt-q33tt)/(3-q11tt-q22tt-q33tt); % Pr[St=1|St-1=1] q11ss=(2-q33tt-q22tt)/(3-q11tt-q22tt-q33tt); % Pr[St=0|St-1=0] prob_T=[q11ss;q11ss;q11ss;q22ss;q22ss;q22ss;q33ss;q33ss;q33ss]; zt_mat=[0;1;1];zt_mat2=[0;0;1]; alfat=alfa3tt*zt_mat2+alfa2tt*zt_mat+alfa1tt*ones(3,1); Ht=[ktt*alfat]; var_L=sigvtt*ones(3,1); filtpr=[]; % Filter for j=1:T; betat=betatt(j,1); yt=yy(j,1)*ones(3,1); fcast=yt-Ht*betat; %form the probability matrix prob_=prob_T.*trmat; prob_T=[prob_(1,1)+prob_(4,1)+prob_(7,1);prob_(2,1)+prob_(5,1)+prob_(8,1);prob_(3,1)+prob_(6,1)+prob_(9,1)]; %Pr[Zt-m,...,Zt-1,Zt|Yt] pr_vl=(1./sqrt(2.*pi.*var_L)).*exp(-0.5*fcast.*fcast./var_L).*prob_T; pr_val=sum(pr_vl); pro_=pr_vl/(pr_val); filtpr=[filtpr prob_T]; %Pr[Zt-m+1,...,Zt-1,Zt|Yt] prob_T=[pro_(1,1);pro_(1,1);pro_(1,1);pro_(2,1);pro_(2,1);pro_(2,1);pro_(3,1);pro_(3,1);pro_(3,1)]; end, % generate ZT q1=filtpr(1,T); % q2=filtpr(2,T); % q3=filtpr(3,T); m=1; %order of the binomial generator Zt=zeros(T,1); pr1=q1/(q1+q2+q3+0.001); pr3=q3/(q1+q2+q3+0.001); u=unifrnd(0,1); Zt(T,1)=le(pr1,u); %binomial generated based on filtered value of Pr[ZT|YT] Zt2(T,1)=ge(u,1-pr3); % generate the rest of Zt for i=2:T-1; istar=T-i; if Zt(istar+1,1)==0; q1=q11tt*filtpr(1,istar); q2=q12tt*filtpr(2,istar); q3=q13tt*filtpr(3,istar); elseif Zt(istar+1,1)==1; if Zt2(istar+1,1)==0; q1=q21tt*filtpr(1,istar); q2=q22tt*filtpr(2,istar); q3=q23tt*filtpr(3,istar); elseif Zt2(istar+1,1)==1; q1=q31tt*filtpr(1,istar); q2=q32tt*filtpr(2,istar); q3=q33tt*filtpr(3,istar); end, end, %generate Zt based on the new probabilities pr1=q1/(q2+q1+q3); pr3=q3/(q2+q1+q3); u=unifrnd(0,1); Zt(istar,1)=le(pr1,u); Zt2(istar,1)=ge(u,1-pr3); end, ZTT=Zt; ZTT2=Zt2; %------------------------------------------------------------------------- % % Function that generates beta(t) (8/30/99) % beta(t) is the unobserved variable variable % Procedure: Extract the parameters and define the state space model % Run the Kalman Filter and save the results % Generate each beta using Cholesky decomposition % beta=[delu etahat] % %------------------------------------------------------------------------- function [betatt]=genxs3d(param,trparamq,trparamp,STT,ZTT,ZTT2); load data48; yy=dat48(206:600,:); T=size(yy,1); % Extract the parameters from the input vector sigett=param(1,1); sigvtt=param(2,1); eta2tt=param(3,1); eta1tt=param(4,1); alfa3tt=param(5,1); alfa2tt=param(6,1); alfa1tt=param(7,1); ktt=param(8,1); q33tt=trparamq(3,3); q22tt=trparamq(2,2); q11tt=trparamq(1,1); q31tt=trparamq(3,1); q32tt=trparamq(3,2); q13tt=trparamq(1,3); q12tt=trparamq(1,2); q23tt=trparamq(2,3); q21tt=trparamq(2,1); p22tt=trparamp(2,2); p11tt=trparamp(1,1); p12tt=trparamp(1,2); p21tt=trparamp(2,1); % calculate the gain sequence vector conditional on generated STT gamat=zeros(T,1); for i=1:T; etat(i,1)=eta2tt*STT(i,1)+eta1tt; %define conditional etas alfat(i,1)=alfa3tt*ZTT2(i,1)+ZTT(i,1)*alfa2tt+alfa1tt; %define condtional alfas end, g=zeros(T,1); g(1,1)=10; gamat(1,1)=1/g(1,1); for i=2:T; if STT(i,1)==STT(i-1,1); %if the state does not change, then lower the gain g(i,1)=g(i-1,1)+1; else, g(i,1)=1; end, gamat(i,1)=1/g(i,1); end, Gt=[1,0;-ktt,1]; % G multiplies the control shock QQ=[sigvtt, 0; 0, sigett]; % covariance matrix for the transition equation RR=[0.01]; % covariance matrix for the measurement equation betattm=zeros(1,T); % storage matrix for beta pttm=[sigvtt]; % storage matrix for the variance %************************************************************************* % Initialize the filter p22ss=(1-p11tt)/(2-p11tt-p22tt+0.0001); p11ss=(1-p22tt)/(2-p11tt-p22tt+0.0001); % Pr[St=1|St-1=1] q33ss=(2-q22tt-q11tt)/(3-q11tt-q22tt-q33tt); q22ss=(2-q33tt-q11tt)/(3-q11tt-q22tt-q33tt); % Pr[St=1|St-1=1] q11ss=(2-q22tt-q33tt)/(3-q11tt-q22tt-q33tt); % Pr[St=0|St-1=0] pre_cf=eta2tt*p22ss+eta1tt; Fss=[0.9]; vecq=[sigvtt;0;0;sigett]; tmp_vr=inv(eye(1)-kron(Fss,Fss))*vecq; pre_vr=tmp_vr(4,1); % Iterate for j=2:T; Ft=[1-gamat(j,1)]; Ht=[ktt*alfat(j,1)]; mumat=[gamat(j,1)*yy(j,2)]; pst_cf=Ft*pre_cf+mumat; %posterior beta pst_vr=Ft*pre_vr*Ft'+RR; %posterior beta variance f_cast=yy(j,1)-Ht*pst_cf; %forecast error ss=Ht*pst_vr*Ht'+(ktt^2)*sigvtt+sigett; %forecast error variance k_gn=pst_vr*Ht'*inv(ss); pre_cf=pst_cf+k_gn*f_cast; pre_vr=(eye(1)-k_gn*Ht)*pst_vr; betattm(:,j)=pre_cf; %store the results pttm=[pttm pre_vr]; end, beta_mat=zeros(T,1); %create the storage matrix for generated betas X_temp=betattm(:,T); %get the last iteration P_temp=pttm(:,T); Xtemp1=X_temp(1,1); ptmp=abs(P_temp(1,1)); GG=sqrt(ptmp); bstar=normrnd(0,1); betat=Xtemp1+GG*bstar; beta_mat(T,1)=betat; for i=1:T-1; istar=T-i; mumat=[gamat(istar,1)*yy(istar,2)]; Ft=[1-gamat(istar,1)]; Ht=[ktt*alfat(istar,1)]; f_cast0=betat-Ft*betattm(:,istar)-mumat; ptemp=pttm(:,istar); ss0=Ft*ptemp*Ft'+RR; k_gn0=ptemp*Ft'/ss0; betattm(:,istar)=betattm(:,istar)+k_gn0*f_cast0; ptemp=(eye(1)-k_gn0*Ft)*ptemp; X_temp=betattm(:,istar); Xtemp1=X_temp(1,1); P_temp=ptemp; ptmp=abs(P_temp(1,1)); GG=sqrt(ptmp); bstar=normrnd(0,1); betat=Xtemp1+GG*bstar; beta_mat(istar,1)=betat; end, betatt=beta_mat; %------------------------------------------------------------------------- % % Function that generates St (10/06/99) % St is the economic state variable indicator % Also generates the conditional probability of each state % Does this generate the transition probability? % %------------------------------------------------------------------------- function [STT]=genst3d(param,trparamq,trparamp,STT,betatt); load data48; yy=dat48(206:600,:); T=size(yy,1); % Extract the parameters from the input vector sigett=param(1,1); sigvtt=param(2,1); eta2tt=param(3,1); eta1tt=param(4,1); alfa3tt=param(5,1); alfa2tt=param(6,1); alfa1tt=param(7,1); ktt=param(8,1); q33tt=trparamq(3,3); q22tt=trparamq(2,2); q11tt=trparamq(1,1); q31tt=trparamq(3,1); q32tt=trparamq(3,2); q13tt=trparamq(1,3); q12tt=trparamq(1,2); q23tt=trparamq(2,3); q21tt=trparamq(2,1); p22tt=trparamp(2,2); p11tt=trparamp(1,1); p12tt=trparamp(1,2); p21tt=trparamp(2,1); STT=zeros(T,1); st_mat=[0;1]; mutt=eta1tt*ones(2,1)+eta2tt*st_mat; var_L=(sigett+(ktt^2)*sigvtt)*ones(2,1); Stm1=0; gt=2; for i=2:T-1; Stp1=STT(i+1,1); yt=yy(i,2)*ones(2,1); % choose the transition probabiity if Stm1==0; probs=[p11tt;p12tt]; else, probs=[p21tt;p22tt]; end, if Stp1==0; prs=[p11tt;p21tt]; else, prs=[p12tt;p22tt]; end, % determine forecast error fcast=yt-mutt; % determine the forecast eror probability pr_v_l=(1./sqrt(2.*pi.*var_L)).*exp(-0.5*fcast.*fcast./var_L); pr=probs.*pr_v_l.*prs; p1=pr(1,1)/sum(pr); p2=pr(2,1)/sum(pr); % now generate u=unifrnd(0,1); St=le(p1,u); STT(i,1)=St; % reset and adjust the gain sequence if St==Stm1; gt=gt+1; else, gt=0; end, Stm1=St; end, STT(T,1)=STT(T-1,1); %------------------------------------------------------------------------- % % Function that generates mu(t) (10/06/99) % mu(t) is the vector that maps the past unobserved variable to the % current unobserved variable. Mu is dependent on the value % of the state variable eta (consequently on the state variable St). % Procedure: This is a function that generates the hyperparameters % etah and etal. eta_t=etah*Stt+etal*(1-Stt) % eta_t forms the vector mu_t = (eta_t 0)' %------------------------------------------------------------------------- function [eta2tt,eta1tt,sigett]=genmu3d(param,trparamq,trparamp,STT,betatt); load data48; yy=dat48(206:600,:); T=size(yy,1); R0_mu=eye(2)/4; % prior for etah and etal t0_mu=[0;0]; d0=0.001; % Extract the parameters from the input vector sigett=param(1,1); sigvtt=param(2,1); eta2tt=param(3,1); eta1tt=param(4,1); alfa3tt=param(5,1); alfa2tt=param(6,1); alfa1tt=param(7,1); ktt=param(8,1); q33tt=trparamq(3,3); q22tt=trparamq(2,2); q11tt=trparamq(1,1); q31tt=trparamq(3,1); q32tt=trparamq(3,2); q13tt=trparamq(1,3); q12tt=trparamq(1,2); q23tt=trparamq(2,3); q21tt=trparamq(2,1); p22tt=trparamp(2,2); p11tt=trparamp(1,1); p12tt=trparamp(1,2); p21tt=trparamp(2,1); Xstar=[ones(T,1),STT(1:T,1)]; Ystar=yy(:,2); V=inv(R0_mu+Xstar'*Xstar); ETA=V*(R0_mu*t0_mu+Xstar'*Ystar); C=chol(V); accept=0; waita=1; while accept==0, etat=ETA+C'*normrnd(0,1,2,1); if etat(2,1) > 0; accept=1; else, waita=waita+1; end, if waita > 1000; etat(2,1)=eta2tt; accept=1; end, end, eta1tt=etat(1,1); eta2tt=etat(2,1); %need to fix this gap=zeros(T,1); gap(1,1)=0; etapr=zeros(T,1); for i=2:T; gap(i,1)=yy(i,2)-(eta1tt+eta2tt*STT(i,1)); end, d=d0+gap'*gap; c=chi2rnd(T); sigett=d/(c+ktt^2); %------------------------------------------------------------------------- % % Function that generates F(t) (7/22/99) % F(t) is the matrix that maps the past unobserved variable to the % current unobserved variable. F is dependent on the value % of the state variable eta (consequently on the state variable % St). In the learning model, F is a fucntion of the time since % the last switch, i.e., a reset of the gain sequence. % Procedure: This is a function that generates the hyperparameters for F % %------------------------------------------------------------------------- function [ktt]=genk3d(param,trparamq,trparamp,STT,ZTT,ZTT2,betatt); load data48; yy=dat48(206:600,:); T=size(yy,1); % Bayesian Priors for calculation of the Hyperparameters R0_k=1; t0_k=0; sigett=param(1,1); sigvtt=param(2,1); eta2tt=param(3,1); eta1tt=param(4,1); alfa3tt=param(5,1); alfa2tt=param(6,1); alfa1tt=param(7,1); ktt=param(8,1); q33tt=trparamq(3,3); q22tt=trparamq(2,2); q11tt=trparamq(1,1); q31tt=trparamq(3,1); q32tt=trparamq(3,2); q13tt=trparamq(1,3); q12tt=trparamq(1,2); q23tt=trparamq(2,3); q21tt=trparamq(2,1); p22tt=trparamp(2,2); p11tt=trparamp(1,1); p12tt=trparamp(1,2); p21tt=trparamp(2,1); z1=ones(T,1); ztt=[z1 ZTT ZTT2]; alfas=[alfa1tt;alfa2tt;alfa3tt]; alphas=ztt*alfas; Xstar=betatt.*alphas; Ystar=yy(:,1); V=inv(R0_k+Xstar'*Xstar); K=V*(R0_k*t0_k+Xstar'*Ystar); C=chol(V); ktt=K+C'*normrnd(0,1,1,1); %------------------------------------------------------------------------- % % Function that generates alfa(t) (10/06/99) % alfa(t) is the parameter that maps the past unobserved variable to the % current observed variable. Alfa is dependent on the value % of the state variable Zt. % Procedure: This is a function that generates the hyperparameters % alfah and alfal. alfa_t=alfah*Ztt+alfal*(1-Ztt) % alfa_t forms the matrix H(sig_t) = (0 alfa_t ; 1 0) % %------------------------------------------------------------------------- function [alfa3tt,alfa2tt,alfa1tt,sigvtt]=genalfas3c(param,trparamq,trparamp,ZTT,ZTT2,betatt); load data48; yy=dat48(206:600,:); T=size(yy,1); % Bayesian Priors for calculation of the Hyperparameters R0_H=eye(3)/9; % prior for alfah and alfal t0_H=[0;0;0]; % Extract the parameters from the input vector sigett=param(1,1); sigvtt=param(2,1); eta2tt=param(3,1); eta1tt=param(4,1); alfa3tt=param(5,1); alfa2tt=param(6,1); alfa1tt=param(7,1); ktt=param(8,1); q33tt=trparamq(3,3); q22tt=trparamq(2,2); q11tt=trparamq(1,1); q31tt=trparamq(3,1); q32tt=trparamq(3,2); q13tt=trparamq(1,3); q12tt=trparamq(1,2); q23tt=trparamq(2,3); q21tt=trparamq(2,1); p22tt=trparamp(2,2); p11tt=trparamp(1,1); p12tt=trparamp(1,2); p21tt=trparamp(2,1); % define ystar and xstar xstar=ktt*ZTT.*betatt; xstar1=ktt*ZTT2.*betatt; Xstar=[betatt*ktt,xstar,xstar1]; Ystar=[yy(:,1)]; V=inv(R0_H+Xstar'*Xstar); ALFA=V*(R0_H*t0_H+Xstar'*Ystar); C=chol(V); accept=0; waita=1; while accept==0, alfat=ALFA+C'*normrnd(0,1,3,1); if alfat(2,1) > 0; if alfat(3,1) > 0; accept=1; else, waita=waita+1; end, else, waita=waita+1; end, if waita > 1000; alfat(1,1)=alfa1tt; alfat(2,1)=alfa2tt; alfat(3,1)=alfa3tt; accept=1; end, end, alfa1tt=alfat(1,1); alfa2tt=alfat(2,1); alfa3tt=alfat(3,1); alfapr=zeros(T,1); for i=1:T; alfapr(i,1)=alfa1tt+alfa2tt*ZTT(i,1)+alfa3tt*ZTT2(i,1); end, %this needs to be fixed d0=0.01; gap=yy(:,1)-ktt*alfapr(i,1).*betatt; d=d0+gap'*gap; c=chi2rnd(T); sigvtt=d/c; %************************************************************************* % GIBBS (10/1/2000) %************************************************************************* % This is the main program for the Gibbs sampling estimation % of the model in Owyang, 1999 for the following state space model: % y_t = H(sig_t) * beta_t + e_t % beta_t = F(S_t) * beta_t-1 + mu_t + G * v_t % y_t = (pi_t delu_t)' % H(sig_t) = (0 alfa_t ; 1 0) % beta_t = (delu_t etahat_t)' % F(S_t) = (0 0 ; gama_t 1-gama_t) % mu_t = (eta_t 0)' % G = (1 -k ; 0 0) % v_t = (epsilon v)' % eta=etah if St=1; eta=etal if St=0 % alfa=alfah if Zt=1; alfa=alfal if Zt=0 %************************************************************************ clear all; % The data will be input by each individual function. % Intialize the Gibbs sampler M=0; %number of iterations before the sampler begins to save results N=110; %number of iterations used to calcutlate the results LL=2; %every LLth value is used capN=N+M; %total number of draws load data48; dat48=dat48(6:600,:); T=size(dat48,1)-200; % Initial Conditions sigett=0.2;sigvtt=0.2; eta2tt=0.3;eta1tt=0.01; alfa3tt=1.5;alfa2tt=0.8;alfa1tt=0.01; ktt=1; q33tt=0.8;q22tt=0.8;q11tt=0.8; q31tt=0.1;q32tt=0.1;q13tt=0.1; q12tt=0.1;q23tt=0.1;q21tt=0.1; p22tt=0.8;p11tt=0.8; p12tt=0.1;p21tt=0.1; STT=zeros(T,1); %the initial vector of states, (S1,S2,...,ST) ZTT=zeros(T,1); %the initial vector of states, (Z1,Z2,...,ZT) ZTT2=zeros(T,1); betatt=zeros(T,1); %the beta is a 2x1--betatt is Tx2 % Storage spaces p22mm=[];p11mm=[]; % Pr[St=1|St-1=1] sigvmm=[]; % variance of the control shock sigemm=[]; % variance of the Phillips curve shock kmm=[]; % slope of the Phillips curve eta2mm=[]; eta1mm=[]; % magnitude of the shock in the low state alfa3mm=[]; % policymaker's preferences in inflation hawk state alfa2mm=[]; alfa1mm=[]; % policymaker's preferences in inflation tolerant state q22mm=[]; q11mm=[]; % Pr[Zt=1|Zt-1=1] q33mm=[]; % Pr[Zt=0|Zt-1=0] STTmm=zeros(T,1); % the state of the economy ZTTmm=zeros(T,1); % the policymaker's preference state ZTT2mm=zeros(T,1); p12mm=[];p21mm=[]; q31mm=[];q32mm=[];q13mm=[];q12mm=[];q23mm=[];q21mm=[]; param=zeros(8,1); trparamp=zeros(3,3); trparamq=zeros(2,2); % Begin iterations for n=1:capN; % set up initial parameter vector disp('sampling iteration #');disp(n); param(1,1)=sigett; param(2,1)=sigvtt; param(3,1)=eta2tt; param(4,1)=eta1tt; param(5,1)=alfa3tt; param(6,1)=alfa2tt; param(7,1)=alfa1tt; param(8,1)=ktt; trparamq(3,3)=q33tt; trparamq(2,2)=q22tt; trparamq(1,1)=q11tt; trparamq(3,1)=q31tt; trparamq(3,2)=q32tt; trparamq(1,3)=q13tt; trparamq(1,2)=q12tt; trparamq(2,3)=q23tt; trparamq(2,1)=q21tt; trparamp(2,2)=p22tt; trparamp(1,1)=p11tt; trparamp(1,2)=p12tt; trparamp(2,1)=p21tt; % begin generating [betatt]=genxs3d(param,trparamq,trparamp,STT,ZTT,ZTT2); [ZTT,ZTT2]=genzt3d(param,trparamq,trparamp,betatt,STT); [STT]=genst3d(param,trparamq,trparamp,STT,betatt); [eta2tt,eta1tt,sigett]=genmu3d(param,trparamq,trparamp,STT,betatt); param(1,1)=sigett; param(3,1)=eta2tt; param(4,1)=eta1tt; [ktt]=genkd(param,trparamq,trparamp,STT,ZTT,ZTT2,betatt); param(8,1)=ktt; [alfa3tt,alfa2tt,alfa1tt,sigvtt]=genalfas3d(param,trparamq,trparamp,ZTT,ZTT2,betatt); param(2,1)=sigvtt; param(5,1)=alfa3tt; param(6,1)=alfa2tt; param(7,1)=alfa1tt; gtemp=[1;2]; % generate the switching matrix ST=STT+ones(size(STT,1),1); ntemp=size(ST,1); mtemp=size(gtemp,1); transmat=zeros(mtemp,mtemp); for n0=2:ntemp; Swt1=ST(n0-1); Swt=ST(n0); transmat(Swt1,Swt)=transmat(Swt1,Swt)+1; end, gtemp=[1;2;3]; % generate the switching matrix mtemp=size(gtemp,1); ZT=ZTT+ZTT2+ones(size(ZTT,1),1); ntemp=size(ZT,1); tranzmat=zeros(mtemp,mtemp); for m0=2:ntemp; Zwt1=ZT(m0-1); Zwt=ZT(m0); tranzmat(Zwt1,Zwt)=tranzmat(Zwt1,Zwt)+1; end, % based on the switching matrix, generate the probabilities p11tt=betarnd(transmat(1,1)+1,transmat(1,2)+0.1); p22tt=betarnd(transmat(2,2)+1,transmat(2,1)+0.1); p12tt=1-p11tt; p21tt=1-p22tt; q11tt=betarnd(tranzmat(1,1)+2,tranzmat(1,2)+tranzmat(1,3)+0.1); q22tt=betarnd(tranzmat(2,2)+2,tranzmat(2,1)+tranzmat(2,3)+0.1); q33tt=betarnd(tranzmat(3,3)+2,tranzmat(3,1)+tranzmat(3,2)+0.1); q21tt=betarnd(tranzmat(2,1)+0.5,tranzmat(2,3)+0.5)*(1-q22tt); q23tt=betarnd(tranzmat(2,3)+0.5,tranzmat(2,1)+0.5)*(1-q22tt); q12tt=betarnd(tranzmat(1,2)+0.5,tranzmat(1,3)+0.5)*(1-q11tt); q13tt=betarnd(tranzmat(1,3)+0.5,tranzmat(1,2)+0.5)*(1-q11tt); q31tt=betarnd(tranzmat(3,1)+0.5,tranzmat(3,2)+0.5)*(1-q33tt); q32tt=betarnd(tranzmat(3,2)+0.5,tranzmat(3,1)+0.5)*(1-q33tt); trparamq(3,3)=q33tt; trparamq(2,2)=q22tt; trparamq(1,1)=q11tt; trparamq(3,1)=q31tt; trparamq(3,2)=q32tt; trparamq(1,3)=q13tt; trparamq(1,2)=q12tt; trparamq(2,3)=q23tt; trparamq(2,1)=q21tt; trparamp(2,2)=p22tt; trparamp(1,1)=p11tt; trparamp(1,2)=p12tt; trparamp(2,1)=p21tt; if n > M; % Store the results STTmm=STTmm+STT; % STTmm=[STTmm Stt]; ZTTmm=ZTTmm+ZTT; % ZTTmm=[ZTTmm Ztt]; ZTT2mm=ZTT2mm+ZTT2; p11mm=[p11mm p11tt]; p12mm=[p12mm p12tt]; p21mm=[p21mm p21tt]; p22mm=[p22mm p22tt]; q31mm=[q31mm q31tt]; q32mm=[q32mm q32tt]; q13mm=[q13mm q13tt]; q12mm=[q12mm q12tt]; q23mm=[q23mm q23tt]; q21mm=[q21mm q21tt]; q22mm=[q22mm q22tt]; sigvmm=[sigvmm sigvtt]; % (3) variance of the control shock sigemm=[sigemm sigett]; % (4) variance of the Phillips curve shock kmm=[kmm ktt]; % (5) slope of the Phillips curve eta2mm=[eta2mm eta2tt]; eta1mm=[eta1mm eta1tt]; % (7) magnitude of the shock in the low state alfa3mm=[alfa3mm alfa3tt]; % (8) policymaker's preferences in inflation hawk state alfa2mm=[alfa2mm alfa2tt]; alfa1mm=[alfa1mm alfa1tt]; % (9) policymaker's preferences in inflation tolerant state q22mm=[q22mm q22tt]; q11mm=[q11mm q11tt]; % (10) Pr[Zt=1|Zt-1=1] q33mm=[q33mm q33tt]; % (11) Pr[Zt=0|Zt-1=0] save ch3prob p22mm p11mm q33mm q22mm q11mm; save ch3params eta2mm eta1mm alfa3mm alfa2mm alfa1mm kmm sigvmm sigemm; save ch3ndprob p12mm p21mm q31mm q32mm q13mm q12mm q23mm q21mm; save ch3states STTmm ZTT2mm ZTTmm; end, end, % end of iterations %************************************************************************** % This program is for the simulation of the model in the Owyang (1999) % The program requires the input parameters taken from the Gibbs % sampling program. This is the two state model with fixed state values. % This program was last revised 10/09/99 % The program generates Irate and Du based on the input parameters %************************************************************************** clear all; z=100; Irates=[]; load ch3prob; load ch3params; load ch3states; load ch3ndprob; T=size(eta2mm,2); etam=mean(eta2mm(100:T)); etal=mean(eta1mm(100:T)); alfah=mean(alfa3mm(100:T)); alfam=mean(alfa2mm(100:T)); alfal=mean(alfa1mm(100:T)); sige=mean(sigemm(100:T)); sigv=mean(sigvmm(100:T)); p00=mean(p11mm(100:T)); p11=mean(p22mm(100:T)); q00=mean(q11mm(100:T)); q11=mean(q22mm(100:T)); q22=mean(q33mm(100:T)); p01=mean(p12mm(100:T)); p10=mean(p21mm(100:T)); q01=mean(q12mm(100:T)); q02=mean(q13mm(100:T)); q10=mean(q21mm(100:T)); q12=mean(q23mm(100:T)); q20=mean(q31mm(100:T)); q21=mean(q32mm(100:T)); k=mean(kmm(100:T)); for i=1:z; disp('montecarlo iteration #'),disp(i); Irate=[];Pibar=[];Du=[]; t=50; % time before the regime change %************************************************************************** % Determine the shocks % eta is the intercept shock to the Phillips curve % epsilon is the noise shock to the Phillips curve % nu is the inflation noise %************************************************************************** shcks=zeros(2,t+T); % the first shock is the Phillips curve shock % the second shock is the inflation noise e1=randn(2,t+T); shcks(1,:)=e1(1,:)*sqrt(sige); shcks(2,:)=e1(2,:)*sqrt(sigv); eta=etam*ones(1,50); % seeds the first t periods with etal Eta=[1]; % this generates the values of eta for i=2:T; prox1=rand(1,1); if Eta(1,i-1)==1; if prox1 <= p11; Eta1=1; else, Eta1=0; end, elseif Eta(1,i-1)==0; if prox1 <= p00; Eta1=0; else, Eta1=1; end, end, Eta=[Eta Eta1]; clear prox1; end for i=1:T; if Eta(1,i)==2; eta1=etah; elseif Eta(1,i)==1; eta1=etam; else, eta1=etal; end, eta=[eta eta1]; end, %************************************************************************** % Determine the policymaker preferences through a Markov process %************************************************************************** Alpha=[];Alpha0=0; for i=1:t+T; prox=rand(1,1); if Alpha0==1; if prox <= q11; Alpha1=1; elseif prox >= 1-q12; Alpha1=2; else, Alpha1=0; end, elseif Alpha0==2; if prox <= q22; Alpha1=2; elseif prox >= 1-q20; Alpha1=0; else, Alpha1=1; end, else, if prox <= q00; Alpha1=0; elseif prox >= 1-q02; Alpha1=2; else, Alpha1=1; end, end, Alpha=[Alpha Alpha1]; Alpha0=Alpha1; clear prox; end, alfa=[]; for i=1:t+T; if Alpha(1,i)==2; alfa1=alfah; elseif Alpha(1,i)==1; alfa1=alfam; else, alfa1=alfal; end, alfa=[alfa alfa1]; end, %************************************************************************* % Begin the simulation %************************************************************************* etahat=[1]; pibar=[]; irate=[]; du=[]; dusq=[]; % for the period before the incidence of the first shock for i=1:t; ppibar=k*etahat(1,i)*alfa(1,i); % set the target level pirate=ppibar+shcks(2,i); % inflation is target plus noise pdu=shcks(1,i)+eta(1,i)-k*shcks(2,i); % current unemployment gap pibar=[pibar ppibar]; irate=[irate pirate]; du=[du pdu]; %dusqt=pdu*pdu'; %dusq=[dusq dusqt]; ssqdu=sum(du); betahat=(1/(i+10))*ssqdu; etahat=[etahat betahat]; clear ppibar;clear pirate;clear pdu; end, clear dusq;dusq=[];j=0;Eta0=0; % for the period after the incidence of the shock for i=1:T; ppibar=k*etahat(1,i+t)*alfa(1,i+t); pirate=ppibar+shcks(2,i+t); pdu=shcks(1,i+t)+eta(1,i+t)-k*shcks(2,i+t); pibar=[pibar ppibar]; irate=[irate pirate]; du=[du pdu]; %dusqt=pdu*pdu'; %dusq=[dusq dusqt]; ssqdu=sum(du); if Eta(1,i) == Eta0; % if eta same, lower the gain j=j+1; betahat=(1/j)*ssqdu; else, j=1; % if eta switches then reset the gain sequence clear du; % and clear the previous history du=[]; end, Eta0=Eta(1,i); etahat=[etahat betahat]; clear ppibar;clear pirate;clear pdu; end, for j=t+1:t+T; itemp=irate(1,j);%dtemp=du(1,j); Irate=[Irate;itemp]; %Du=[Du;dtemp]; end, Irates=[Irates Irate]; %clear irate;clear pibar;clear du;clear shcks;clear alpha; end, Irate1=Irates(:,1:200); Irate2=Irates(:,201:400); Irate3=Irates(:,401:600); Irate4=Irates(:,601:800); Irate5=Irates(:,801:1000); save run2_1c Irate1 -ascii; save run2_2c Irate2 -ascii; save run2_3c Irate3 -ascii; save run2_4c Irate4 -ascii; save run2_5c Irate5 -ascii;