clear all

%read in data from Jorgensen KLEM

% %A=dlmread('G:\Research\AllenT\NK Philips Curve\Data\35sector.klem.tab', '\t',0,0);
% B=zeros(46,104,35);
% count=0;
% for j=1:35
%     count=count+1;
%     B(:,:,j)=A( ( (j-1)*46) + 1:j * 46, : );
% end

A=dlmread('G:\Research\AllenT\NK Philips Curve\Data\35klem96.dat', '\t',0,0);
B=zeros(39,13,35);
index=1;
for j=1:35
    B(:,1:5,j) = A(index:index+38, 1:5);
    index=index+39;
    B(:,6:9,j) = A(index:index+38,3:6);
    index=index+39;
    B(:,10:13,j) =  A(index:index+38,3:6);
    index=index+39;
end
    
%X(inflation, labour share, log(labour share), mean(log(labour share)),
%  rmc, rmc adj. for price staggering)
%get prices and labour share for each sector

X=zeros(38,5,35);
X(:,1,:)= log(B(2:39,4,:) ./ (B(1:38,4,:))) ;
X(:,2,:)=B(2:39,12,:) ./ (B(2:39,3,:) .* B(2:39,4,:));

%compute the log deviation of labour share from the mean to get 
%a measure of real marginal costs
%NOTE: I calculate the mean from each sector NOT the whole sample as
%I have no guidance on whether to use the whole or sector

X(:,3,:)= log(X(:,2,:));
a_j = 1-reshape(mean(X(:,2,:)),35,1);
X(:,4,:)= reshape((reshape(log(mean(X(:,2,:))),35,1)*ones(1,38))',38,1,35);
X(:,5,:)= X(:,3,:)-X(:,4,:);
markup= .1;
eta= (1+markup)/(2+markup);
h_j= (eta.*a_j)./(1-a_j);
X(:,6,:)= repmat(1./(1+h_j),1,38)'.*reshape(X(:,5,:),38,35);

%dynare stuff
inf=X(:,1,1);
mc=X(:,6,1);

%initialise some parameters
f=size(X,3);
beta=0.99;
periods=size(X,1);
k=5;
xmatrix=zeros(1,8);

%calculate the weights of each sector, assuming the economy is made up of
%the sectors represented in the data, and the mean share over time

weights = reshape(mean((B(2:39,3,:) .* B(2:39,4,:)),1),35,1) ./ sum(reshape(mean((B(2:39,3,:) .* B(2:39,4,:)),1),35,1));

%allow for different measures of RMC
%set adj to 6 for adj RMC, 5 for RMC, and 2 for plain ol' labour share

adj=6;

%estimate sectoral parameters


ML_matrix=zeros(f,k);
ML_se = zeros(f,k);
for j=1:f
    
    y=X(2:periods,1,j);
    x=[X(1:periods-1,1,j) X(2:periods,adj,j) X(1:periods-1,adj,j)];
    y1=X(3:periods,adj,j);
    x1=[X(2:periods-1,adj,j) X(1:periods-2,adj,j)];
    
    [xhat,se] = lscov(x,y);
    [xhat2,se2] = lscov(x1,y1);
    
    
    ML_matrix(j,1:3)=xhat';
    ML_matrix(j,4:5)=xhat2';
    
    ML_se(j,1:3)=se';
    ML_se(j,4:5)=se2';
 
end
    
%create an aggregate inflation series by aggregating over the sectors

weighted_X = zeros(38,6,35);

for j=1:f
    
    weighted_X(:,:,j) = weights(j) .* X(:,:,j);

end

agg_X = sum(weighted_X,3);
  
%ML estimate of aggregate inflation
theta_hats=zeros(3,5);
sigma_hats=zeros(3,5);

[theta_hats(1,1:3),sigma_hats(1,1:3)]= lscov([agg_X(1:37,1) agg_X(2:38,adj) agg_X(1:37,adj)],agg_X(2:38,1));
[theta_hats(1,4:5),sigma_hats(1,4:5)]= lscov([agg_X(2:periods-1,adj) agg_X(1:periods-2,adj)],agg_X(3:periods,adj));

%MG Estimator

theta_hats(2,:)=weights'*ML_matrix;

%generate the variance/covariance matrix for the MG estimator
%do it for both equations 

    Z=ML_matrix-repmat(theta_hats(2,:),35,1);
    Z2=Z(:,4:5);
    Z(:,4:5)=[];
    
    MGvcv=zeros(3,3);
    for j=1:f
        vcv_j=Z(j,:)'*Z(j,:);
        MGvcv=MGvcv+(1/(f*(f-1)))*vcv_j;
    end
        
    MGvcv2=zeros(2,2);
    for j=1:f
        vcv_j2=Z2(j,:)'*Z2(j,:);
        MGvcv2=MGvcv2+(1/(f*(f-1)))*vcv_j2;
    end
    
sigma_hats(2,1:3)=diag(MGvcv);
sigma_hats(2,4:5)=diag(MGvcv2);
    
%create Z array which allows Z'Z for later, z dimension is firms
%remembering that we are really doing two estimations

    Z=zeros(periods-1,3,f);
    Z(:,1,:)=X(1:periods-1,1,:);
    Z(:,2,:)=X(2:periods,adj,:);
    Z(:,3,:)=X(1:periods-1,adj,:);
    
    Z2=zeros(periods-2,2,f);
    Z2(:,1,:)=X(2:periods-1,adj,:);
    Z2(:,2,:)=X(1:periods-2,adj,:);  

%create estimate of sigma for RC estimation
    sigma=zeros(f,1);
    sum_est=zeros(3,3,f);
    for j=1:f
        sigma(j)=(1/(periods-3))*X(2:periods,1,j)'*(eye(periods-1)-Z(:,:,j)*inv(Z(:,:,j)'*Z(:,:,j))*Z(:,:,j)')*X(2:periods,1,j);      
        sum_est(:,:,j)=sigma(j)*inv(Z(:,:,j)'*Z(:,:,j));
    end
    
    sigma2=zeros(f,1);
    sum_est2=zeros(2,2,f);
    for j=1:f
        sigma2(j)=(1/(periods-2))*X(3:periods,adj,j)'*(eye(periods-2)-Z2(:,:,j)*inv(Z2(:,:,j)'*Z2(:,:,j))*Z2(:,:,j)')*X(3:periods,adj,j);      
        sum_est2(:,:,j)=sigma2(j)*inv(Z2(:,:,j)'*Z2(:,:,j));
    end
    

%generate the RC weights
    delta=f*MGvcv - mean(sum_est,3);
    W=zeros(3,3,f);
    %sum(,3)
    var_j=zeros(3,3,35);
    for h=1:f
         var_j(:,:,h)=inv(delta+sum_est(:,:,h));
    end
    term=sum(var_j,3); 
    term=inv(term);
    for j=1:f  
        W(:,:,j)=term*var_j(:,:,j);
    end
    
    delta2=f*MGvcv2 - mean(sum_est2,3);
    W2=zeros(2,2,f);
    var_j2=zeros(2,2,35);
    for h=1:f
         var_j2(:,:,h)=inv(delta2+sum_est2(:,:,h));
    end
    term2=sum(var_j2,3); 
    term2=inv(term2);
    for j=1:f  
        W2(:,:,j)=term2*var_j2(:,:,j);
    end

%generate the RC esimator 
    OLS=reshape(ML_matrix(:,1:3)',1,3,f);
    weighted_est=zeros(3,1,f);
    for j=1:f
        weighted_est(:,:,j)=W(:,:,j)*OLS(:,:,j)';
    end;
    theta_hats(3,1:3)=sum(weighted_est,3);
    sigma_hats(3,1:3)=diag(term);
    
    OLS2=reshape(ML_matrix(:,4:5)',1,2,f);
    weighted_est=zeros(2,1,f);
    for j=1:f
        weighted_est(:,:,j)=W2(:,:,j)*OLS2(:,:,j)';
    end;
    theta_hats(3,4:5)=sum(weighted_est,3);
    sigma_hats(3,4:5)=diag(term2);
    
%create the expected inflation series, i.e E(t-1) Pi_t

    exp_inf = zeros(periods,f);
    for j=1:f
        exp_inf(3:periods,j) = ML_matrix(j,1).*X(2:periods-1,1) + ...
            ML_matrix(j,2).*(ML_matrix(j,4).*X(2:periods-1,adj) + ...
            ML_matrix(j,5).*X(1:periods-2,adj)) + ...
            ML_matrix(j,3).*X(2:periods-1,adj);
    end
    
%estimate the reduced form equation for each sector
ML_matrix2=zeros(f,3);
for j=1:f
   
    y=X(2:periods-1,1,j);
    x=[X(1:periods-2,1,j) exp_inf(3:periods,j) X(2:periods-1,adj,j)];
    
    [xhat,se] = lscov(x,y);
    ML_matrix2(j,:) =xhat'; 
end    

%generate the aggregate, this shouldn't be different from the other
%aggregation, IT ISNT, THEREFORE THIS PART IS REDUNDANT but I'll leave it
%as a check
weighted_X = zeros(36,3,35);

for j=1:f
    
    weighted_X(:,:,j) = weights(j) * [X(1:periods-2,1,j) exp_inf(3:periods,j) X(2:periods-1,adj,j)];

end

agg_X2 = sum(weighted_X,3);

%estimate ML on the aggregate
periods=36;
[xhat,se] = lscov([agg_X2(1:periods-2,1) agg_X2(3:periods,2) agg_X2(2:periods-1,3)],agg_X2(2:periods-1,1));

theta_hats2=zeros(3,3);
sigma_hats2=zeros(3,3);

theta_hats2(1,:) = xhat;
sigma_hats2(1,:) = se;
%estimate the MG for the aggregate

theta_hats2(2,:)=weights'*ML_matrix2;

%generate the variance/covariance matrix for the MG estimator
%do it for both equations 

    Z=ML_matrix2-repmat(theta_hats2(2,:),35,1);
    
    MGvcv=zeros(3,3);
    for j=1:f
        vcv_j=Z(j,:)'*Z(j,:);
        MGvcv=MGvcv+(1/(f*(f-1)))*vcv_j;
    end
        
sigma_hats2(2,1:3)=diag(MGvcv);

%create Z array which allows Z'Z for later, z dimension is firms
%remembering that we are really doing two estimations

    Z=zeros(periods-2,3,f);
    Z(:,1,:)=X(1:periods-2,1,:);
    Z(:,2,:)=exp_inf(3:periods,:);
    Z(:,3,:)=X(2:periods-1,adj,:);

%create estimate of sigma for RC estimation
    periods=size(Z,1);
    sigma=zeros(f,1);
    sum_est=zeros(3,3,f);
    for j=1:f
        sigma(j)=(1/(periods-3))*X(2:periods+1,1,j)'*(eye(periods)-Z(:,:,j)*inv(Z(:,:,j)'*Z(:,:,j))*Z(:,:,j)')*X(2:periods+1,1,j);      
        sum_est(:,:,j)=sigma(j)*inv(Z(:,:,j)'*Z(:,:,j));
    end

%generate the RC weights

    delta=f*MGvcv - mean(sum_est,3);
    W=zeros(3,3,f);
    %sum(,3)
    var_j=zeros(3,3,35);
    for h=1:f
         var_j(:,:,h)=inv(delta+sum_est(:,:,h));
    end
    term=sum(var_j,3); 
    term=inv(term);
    for j=1:f  
        W(:,:,j)=term*var_j(:,:,j);
    end
    
    delta2=f*MGvcv2 - mean(sum_est2,3);
    W2=zeros(2,2,f);
    var_j2=zeros(2,2,35);
    for h=1:f
         var_j2(:,:,h)=inv(delta2+sum_est2(:,:,h));
    end
    term2=sum(var_j2,3); 
    term2=inv(term2);
    for j=1:f  
        W2(:,:,j)=term2*var_j2(:,:,j);
    end

%generating the RC estimator

%compare estimators across simulations

theta_hats;
sigma_hats;

theta_hats2;
sigma_hats2;