% Set up data and calendars for monthly data in pcomp .. convert data to quarterly before using
% 1/23/2015, mww

if load_data == 1;

% ----------- Features of Data Set ---------
ns = 26;
miss_code = 1.0e+32;

% ----------- Sample Period, Calendars and so forth
[dnobs_m,calvec_m,calds_m] = calendar_make([1959 1],[2015 3],12);
[dnobs_q,calvec_q,calds_q] = calendar_make([1959 1],[2015 1],4);



% --------------- Read In Quarterly Data for Nominal PCE ---------------- 
xlsname = '/Users/mwatson/Dropbox/p_comp/ddisk/data/nom_pce_m59.xlsx';
% Read Data 
ndesc=1;
ncodes=1;
sheet='Sheet1';
[namevec,descmat,tcodemat,datevec,datamat_m] = readxls(xlsname,sheet,ns,dnobs_m,ndesc,ncodes);
tcodemat = tcodemat(1,:);
labelvec=descmat(:,1);
% Convert Namestrings to upper case 
namevec = upper(namevec);
% Eliminate any leading or trailing blanks 
namevec=strtrim(namevec);
labelvec=strtrim(labelvec);
% Replace missing values with NaN
isel = datamat_m == miss_code;
datamat_m(isel) = NaN;
datamat = mtoq(datamat_m,calds_m,calds_q,0);
pce_nom = datamat;
namevec_nom = namevec;
labelvec_nom = labelvec;
tcodemat_nom = tcodemat;

% ------------- Read in Price Deflators and compute inflation Rates ---------
%xlsname = '/Users/mwatson/Dropbox/p_comp/data/p_pce_work1.xlsx';
xlsname = '/Users/mwatson/Dropbox/p_comp/data/p_pce_m59.xlsx';
% Read Data 
ndesc=2;
ncodes=2;
sheet='Sheet1';
[namevec,descmat,codemat,datevec,datamat_m] = readxls(xlsname,sheet,ns,dnobs_m,ndesc,ncodes);
tcodemat = codemat(1,:);
ordermat = codemat(2,:)';
labelvec=descmat(:,1);
labelvec_short = descmat(:,2);
% Convert Namestrings to upper case 
namevec = upper(namevec);
% Eliminate any leading or trailing blanks 
namevec=strtrim(namevec);
labelvec=strtrim(labelvec);
labelvec_short=strtrim(labelvec_short);
% Replace missing values with NaN
isel = datamat_m == miss_code;
datamat_m(isel) = NaN;
datamat = mtoq(datamat_m,calds_m,calds_q,0);
pce_p = datamat;
namevec_p = namevec;
labelvec_p = labelvec;
labelvec_short_p = labelvec_short;
tcodemat_p = tcodemat;

% Use Quarterly Calendar
dnobs = dnobs_q;
calvec = calvec_q;
calds = calds_q;

% Compute Total and Values for Included Series
isel = tcodemat_p==2;
p_agg = pce_p(:,isel==1);         % Aggregrate Price Indenx  
isel = tcodemat==2;
nom_agg = pce_nom(:,isel==1);     % Aggregate Nominal PCE
isel = tcodemat_p==5;
p_agg_xfe = pce_p(:,isel==1);     % Aggregrate Price Index XFE  
isel = tcodemat==5;
nom_agg_xfe = pce_nom(:,isel==1);  % Aggregate Nominal PCE XFE

% Compute Energy Services Component of Household Expenditures
isel = tcodemat==3;
nom_eg = pce_nom(:,isel==1);   % Gasoline and other energy goods
isel = tcodemat==4;
nom_hu = pce_nom(:,isel==1);    % Housing and utilities
isel = tcodemat==6;
nom_egs = pce_nom(:,isel==1);   % Energy Goods and Services
isel = tcodemat_p==3;
p_eg = pce_p(:,isel==1);       % Gasoline and other energy goods
isel = tcodemat_p==4;
p_hu = pce_p(:,isel==1);        % Housing and utilities
isel = tcodemat_p==6;
p_egs = pce_p(:,isel==1);       % Energy Goods and Services

nom_hu_es = nom_egs - nom_eg;   % Energy services component of housing and utilities
nom_hu_xes = nom_hu - nom_hu_es; % Housing and Utilities excluding energy services

% Compute some inflation rates
infl_eg = NaN*zeros(dnobs,1);
infl_eg(2:end) = log(p_eg(2:end)./p_eg(1:end-1));
infl_egs = NaN*zeros(dnobs,1);
infl_egs(2:end) = log(p_egs(2:end)./p_egs(1:end-1));
s_eg_egs = nom_eg./nom_egs;   % Share of energy goods in energy goods and services 
infl_hu_es = (infl_egs - (s_eg_egs.*infl_eg))./(1-s_eg_egs);
p_hu_es = NaN*zeros(dnobs,1);
p_hu_es(1) = 100;
for i = 2:dnobs;
    p_hu_es(i) = p_hu_es(i-1)*exp(infl_hu_es(i));
end;
infl_hu = NaN*zeros(dnobs,1);
infl_hu(2:end) = log(p_hu(2:end)./p_hu(1:end-1));
s_es_hu = nom_hu_es./nom_hu;  % Share of energy in housing and utilities
infl_hu_xes = (infl_hu - (s_es_hu.*infl_hu_es))./(1-s_es_hu);
p_hu_xes = NaN*ones(dnobs,1);
p_hu_xes(1) = 100;
for i = 2:dnobs;
    p_hu_xes(i) = p_hu_xes(i-1)*exp(infl_hu_xes(i));
end;

% % Check Calculations
% infl_hu1 = NaN*ones(dnobs,1);
% infl_hu1(2:end) = s_es_hu(2:end).*(log(p_hu_es(2:end)./p_hu_es(1:end-1))) + (1-s_es_hu(2:end)).*(log(p_hu_xes(2:end)./p_hu_xes(1:end-1)));
% infl_egs1 = NaN*ones(dnobs,1);
% infl_egs1(2:end) = s_eg_egs(2:end).*(log(p_eg(2:end)./p_eg(1:end-1))) + (1-s_eg_egs(2:end)).*(log(p_hu_es(2:end)./p_hu_es(1:end-1)));

% Add New Series to various matrices
pce_nom = [pce_nom nom_hu_es nom_hu_xes];
tcodemat = [tcodemat 1 1];
namevec = [namevec;'nom_hu_es';'nom_hu_xes'];
labelvec = [labelvec;'Energy Services in Housing and Utilities';'Housing and Utilities excl Energy Services'];

pce_p = [pce_p [p_hu_es p_hu_xes]];
tcodemat_p = [tcodemat_p 1 1];
ordermat_p = [ordermat;9.2;9.1];
namevec_p = [namevec_p;'p_hu_es';'p_hu_xes'];
labelvec_p = [labelvec_p;'Energy Services in Housing and Utilities';'Housing and Utilities excl Energy Services'];
labelvec_short_p = [labelvec_short_p;'Housing Energy Serv.';'Housing (excl. Energy)'];


isel = (tcodemat_p==1)+(tcodemat_p==3);
p_disagg = pce_p(:,isel==1);
ordervec = ordermat_p(isel'==1);
nom_disagg = pce_nom(:,isel==1);
namevec_disagg = namevec(isel==1);
labelvec_disagg = labelvec_p(isel==1);
labelvec_short_disagg = labelvec_short_p(isel==1);
n_incl = size(p_disagg,2);

% ... Data Check ... 
 % Add nominal disaggegrated PCE and compare to aggregate
 tmp = sum(nom_disagg,2);
 dif = (tmp - nom_agg)./nom_agg; 
 if max(abs(dif)) > .0001;
     error('nominal dissaggregates do not add to nominal aggregate');
 end;
 
 % ... Compute Inflation for Aggregate and Dissagregate PCE
 dp_agg = NaN*ones(dnobs,1);
 dp_agg_xfe = NaN*ones(dnobs,1);
 dp_disagg = NaN*ones(dnobs,n_incl);
 dp_agg(2:dnobs,:) = 400*log(p_agg(2:end,:)./p_agg(1:end-1,:));  
 dp_agg_xfe(2:dnobs,:) = 400*log(p_agg_xfe(2:end,:)./p_agg_xfe(1:end-1,:));  
 dp_disagg(2:dnobs,:) = 400*log(p_disagg(2:end,:)./p_disagg(1:end-1,:));
 
 
%  % Plot Series
%  plot(calvec,dp_agg);
%  title('Aggregate Inflation');
%  waitforbuttonpress;
%  
%  for i = 1:n_incl;
%      plot(calvec,dp_disagg(:,i));
%      str = labelvec_disagg(i);
%      title(str);
%      waitforbuttonpress;
%  end;

 
 %  Compute Shares of nominal PCE
 tmp = sum(nom_disagg,2);
 share_disagg = nom_disagg./repmat(tmp,1,n_incl);
 share_avg = NaN*zeros(dnobs,n_incl);
 share_avg(2:end,:) = (nom_disagg(2:end,:)+nom_disagg(1:end-1,:))./(repmat(tmp(2:end),1,n_incl)+repmat(tmp(1:end-1),1,n_incl));
 
%  % Plot Shares
%  for i = 1:n_incl;
%      plot(calvec,[share_disagg(:,i) share_avg(:,i)]);
%      str = labelvec_disagg(i);
%      title(str);
%      waitforbuttonpress;
%  end;
 

 % Compute Approximations to Aggregate inflation
 tmp = dp_disagg.*share_avg;
 dp_agg_approx1 = sum(tmp,2);
 weight_avg = mean(share_disagg,1)';
 dp_agg_approx2 = dp_disagg*weight_avg;
 dif1 = dp_agg - dp_agg_approx1;
 dif2 = dp_agg - dp_agg_approx2;
 

%  plot(calvec,[dp_agg dp_agg_approx1 dif1]);
%  waitforbuttonpress;
%  plot(calvec,[dp_agg dp_agg_approx2 dif2]);


% Reorder components by ordervec;
%[tmp,i_order] = sort(weight_avg,'descend');
[tmp,i_order] = sort(ordervec,'ascend');
nom_disagg = nom_disagg(:,i_order);
dp_disagg = dp_disagg(:,i_order);
share_avg = share_avg(:,i_order);
weight_avg = weight_avg(i_order);
labelvec_disagg = labelvec_disagg(i_order);
labelvec_short_disagg = labelvec_short_disagg(i_order);
namevec_disagg = namevec_disagg(i_order);

% Compute dp_agg_xe:  aggregate inflation excluding energy
s_egs_agg = nom_egs./nom_agg;   % Share of energy goods and services in total
dp_egs = 400*infl_egs;          % Energy goods and services Inflation
dp_agg_xe = (dp_agg - s_egs_agg.*dp_egs)./(1-s_egs_agg);

% Compute a new dissagregated version of inflation using cor_xfe, food, and
% energy
str='DFXARG';    % food
str=upper(str);
j = colnumber(str,namevec_disagg);
dp_food = dp_disagg(:,j);
nom_food = nom_disagg(:,j);
isel = tcodemat==5;
nom_agg_xfe = pce_nom(:,isel==1);     % Nominal PCE 
nom_disagg_3comp = [nom_agg_xfe nom_food nom_egs];
tmp = sum(nom_disagg_3comp,2);
share_disagg_3comp = nom_disagg_3comp./repmat(tmp,1,3);
share_avg_3comp = NaN*zeros(dnobs,3);
share_avg_3comp(2:end,:) = (nom_disagg_3comp(2:end,:)+nom_disagg_3comp(1:end-1,:))./(repmat(tmp(2:end),1,3)+repmat(tmp(1:end-1),1,3));
weight_avg_3comp = mean(share_disagg_3comp,1)';
dp_disagg_3comp = [dp_agg_xfe dp_food dp_egs];
labelvec_disagg_3comp = ['CoreXFE';'Food   ';'Energy '];
namevec_disagg_3comp = ['CoreXFE';'Food   ';'Energy '];
n_incl_3comp = 3;

 % Save Variable Series 
  slist = {...
           'nom_agg' ... 
           'nom_disagg' ...   
           'dp_agg' ...
           'dp_agg_xfe' ...
           'dp_agg_xe' ...
           'dp_disagg' ...
           'dp_agg_approx1' ...
           'dp_agg_approx2' ...
           'share_avg' ...
           'weight_avg' ...
           'calvec' ...
           'calds' ...
           'dnobs' ...
           'n_incl' ...
           'labelvec_disagg' ...
           'labelvec_short_disagg' ...
           'namevec_disagg' ...
           'dp_disagg_3comp' ...
           'share_avg_3comp' ...
           'weight_avg_3comp' ...
           'labelvec_disagg_3comp' ...
           'namevec_disagg_3comp' ...
           'n_incl_3comp' ...
          }';
  str_tmp = [matdir 'slist'];
  save(str_tmp,'slist');
  for iseries = 1:size(slist,1);
     ustr = char(slist(iseries));
     str_tmp = [matdir ustr];
     str_tmp1 = 'save(str_tmp,''';
     str_tmp2 = ''');';
     eval([str_tmp1 ustr str_tmp2]);
  end;

  
end;  %Load Data end;
 
  
% Load Variables and Give Standard Names 
str_tmp = [matdir 'slist'];
load(str_tmp); 
for iseries = 1:size(slist,1);
     ustr = char(slist(iseries));
     str_tmp = [matdir ustr];
     load(str_tmp);
end; 