%create_structure.m creates the matrices that constitute the structure of
%the semi-structural model
%The system in state-space form is 
%xi_t+1 = F xi_t + PSI eps_t
%o_t = A' + H'xi_t + mu_t
%E(PSI eps'ep'PSI') =PSI*PSI'= Q;
%E(mu,mu') = R;
%where xi_t is a vector of unobservable states, y_t is a vector of observables, and 
% x_t is a vector of exogenous/predetermined variables (in the present model, x_t contains lagged values of Y_t), and e is an intercept. 
%Replication materials for: ``What Do Long Data Tell Us About the Permanent Component of Inflation?''
%© Stephanie Schmitt-Grohé and Martín Uribe, May 2022. 

%Estimation sample
year0 = 1900;
year1 = 2022;
disp([year0 year1])

ny = 3; %# of latent endogenous variables 
L = 1; %# of lags 
nlags = max(L, 2); 
nobs = 3; %number of observable variables
nshocks = 5;

RHO= eye(nshocks);  %matrix of autoregressive coefficients of exogenous bloc
%Distinguish shocks with high and low serial correlation
aux = [1 2 3 5]';
aux_RHOlow = sub2ind(size(RHO),aux,aux);
aux = 4;
aux_RHOhigh = sub2ind(size(RHO),aux,aux);

P = eye(nshocks);
aux_P = find(P);

B = zeros(ny,ny*nlags); %matrix of autoregressive coefficients of endogenous bloc
aux_Bown = find(eye(ny));
aux_Bother = setdiff((1:numel(B(:,1:ny)))',aux_Bown);

%Columns: dXm zm dXn zn dXr
%Rows yhat paihat ihat
C = zeros(ny,nshocks);
%normalize C(3,2)=C(1,4)=1
aux = [3 2; 1 4];
aux_Cones = sub2ind(size(C),aux(:,1),aux(:,2));
C(aux_Cones) = 1;
%zero  impact effect of zm on y or pai. 
aux = [1 2; 2 2];
aux_Czero = sub2ind(size(C),aux(:,1),aux(:,2));
C(aux_Czero) = 0;
%Prior of zero impact effect of Xm on pai  and i. Requires imposing a mean of -1 on elements (2,1) and (3,1)
aux = [2 1; 3 1];
aux_Cminus1 = sub2ind(size(C),aux(:,1),aux(:,2));
%The remaining elements of C have a normal (0,1) prior distribution
aux_Cnormal = setdiff((1:numel(C))',[aux_Cones;aux_Czero;aux_Cminus1]);

Ap  = zeros(nobs,1);
aux_Ap = 1;

R = zeros(nobs);
aux = [(1:nobs)'  (1:nobs)' ];
aux_R = sub2ind(size(R),aux(:,1),aux(:,2));

Hp = zeros(nobs,ny*nlags+nshocks);
%Dy
Hp(1,1) = 1;
Hp(1,ny+1) = -1;
Hp(1,ny*nlags+3) = 1;

%Dpai
Hp(2,2) = 1;
Hp(2,ny+2) = -1;
Hp(2,ny*nlags+1) = 1;

%Di
Hp(3,3) = 1;
Hp(3,ny+3) = -1;
Hp(3,ny*nlags+1) = 1;
Hp(3,ny*nlags+5) = 1;

%The single estimated element of Hp is the one controlling the permanent effect of permanent real-interest-rate shocks on output
aux = [1 ny*nlags+5];
aux_Hp = sub2ind(size(Hp),aux(:,1),aux(:,2));

%Create vector  of indices mapping elements of the vector collecting all of the estimated parameters to the estimated parameters of the specific matrices of the model. 
%For example iCnormal is defined as follows: 
% if p is a vector containing values for all estimated parameters, then the value of  the parameters  C(aux_Cnormal) is p(iCnormal). 
pname = {'RHOlow';'RHOhigh';'P';'Bown';'Bother';'Cnormal';'Cminus1';'Ap';'R';'Hp'};
n2=0;
np = 0; %number of estimated parameters
for i=1:length(pname)
%number of estimated parameters 
eval(['np = np + numel(aux_' pname{i} ');'])
%number of parameters in each estimated subset of parameters: for example, nCnormal  is the number of elements of aux_Cnormal
eval(['n' pname{i} ' = numel(aux_' pname{i} ');']);
%Create iCnormal, etc
n1=n2+1;eval(['n2=n2+n'  pname{i} ';']); 
eval(['i' pname{i} ' =transpose(n1:n2);'])
end%for i=1:length(pname)

%prior probability distributions
bppd = 1; %beta
gppd = 2; %gamma
nppd = 3; %normal
uppd = 4; %uniform 

%parameters of prior distributions
%mean std.dev. min max
ppd_RHOlow  = (iRHOlow*0+1)*[0.3 0.2 0 1 bppd];

ppd_RHOhigh = (iRHOhigh*0+1)*[0.5 0.2 0 1 bppd];

ppd_P =(iP*0+1)* [1 1 0 inf gppd];
ppd_Bown =(iBown*0+1)* [0.5 0.3 -inf inf nppd];
ppd_Bother = (iBother*0+1)*[0 0.3 -inf inf nppd];
ppd_Cnormal =(iCnormal*0+1)* [0 1 -inf inf nppd];
ppd_Cminus1 =(iCminus1*0+1)* [-1 1 -inf inf nppd]; 


%read data
if ~exist('varY') | ~exist('year') | isempty(year)
[Y,year,y,pai,irate,dy,dpai,dirate] = trim_data(year0,year1);
Y = [dy dpai dirate];
T = size(Y,1);
meanY = mean(Y)';
stdY = std(Y)';
varY = var(Y)';
end%if ~exist('varY') | ~exist('year') | isempty(year)
ppd_Ap = [meanY(aux_Ap) stdY(aux_Ap)/sqrt(T)  -inf*abs(meanY(aux_Ap)) inf*abs(meanY(aux_Ap)) meanY(aux_Ap)*0+nppd];
ppd_R = [0*varY  0.1*varY 0*varY  0.1*varY 0*varY+uppd];
ppd_R(:,1) = (ppd_R(:,3)+ppd_R(:,4))/2;
ppd_R(:,2) = (ppd_R(:,4)-ppd_R(:,3))/sqrt(12);
ppd_Hp = [0 5 -inf inf nppd];

%collect prior probability distributions in matrix PPD
PPD = [];
for i=1:length(pname)
eval(['PPD(i' pname{i} ',:) = ppd_' pname{i} ';'])
end%for i=1:length(pname)

%Consolidate names of parameter matrices
aux_C = [aux_Cnormal; aux_Cminus1];
iC = [iCnormal; iCminus1];
[aux_C,aux] = sort(aux_C);
iC = iC(aux);

aux_B = [aux_Bown;aux_Bother];
iB = [iBown;iBother];
[aux_B,aux] = sort(aux_B);
iB = iB(aux);

aux_RHO = [aux_RHOlow;aux_RHOhigh];
iRHO = [iRHOlow;iRHOhigh];
[aux_RHO,aux] = sort(aux_RHO);
iRHO = iRHO(aux);

%Pname is just info: it gives the matrix of parameters to which each element of p belongs
aux = 0; 
for i=1:length(pname);
eval(['J=length(i' pname{i} ');'])
for j=1:J
aux = aux+1;
Pname{aux,1}=pname{i};
end%for j=1:J
end%for i=1:length(pname);


X = 0*Y(:,1)+1; %intercept 

