%SETUPIO   Start with BEA's 2007 use of commodities by industries table 
%          (before redefinitions at producers' prices) and make
%          the following adjustments:
%
%             1. Add used and noncomparable imports to commodity 31G
%             2. Recategorize certain intermediate inputs (e.g., 
%                advertising) as intangible investments
%             3. Replace investment (fixed plus inventories) with
%                two capital flow tables--tangibles and intangibles
%             4. Construct parameter inputs for multi-sector model.
%                       
%          Note: the capital flow table for intangibles differs 
%                conceptually from the BEAs.  R&D is mapped to the
%                commodity of the funder (e.g., pharmacueticals) 
%                rather than to R&D services in NAICS 541.
%    

%          Ellen McGrattan, 5-30-15
%          Revised, ERM, 6-20-16

load UseBPro07.dat 
load detail07.dat

x          = UseBPro07;
y          = detail07;

Iused      = x(16,1:15);
Inonc      = x(17,1:15);
Fused      = x(16,18:23);
Fnonc      = x(17,18:23);

Intbea     = x(1:15,1:15);
WL         = x(20,1:15);
TX         = x(21,1:15);
RKbea      = x(22,1:15);
VAbea      = x(24,1:15);
IObea      = x(25,1:15);

Cbea       = x(1:15,18);
Xfbea      = x(1:15,19);
Xinv       = x(1:15,20);
EX         = x(1:15,21);
IM         = -x(1:15,22);
Gbea       = x(1:15,23);
FUbea      = x(1:15,26);
CObea      = x(1:15,27);

GDPbea     = x(24,26);
GObea      = x(25,27);

%
% Combine industries in detailed IO table
%

y(:,1)     = sumnan(detail07(:,1:13)')';        % 11
y(:,2)     = sumnan(detail07(:,14:21)')';       % 21
y(:,3)     = sumnan(detail07(:,22:24)')';       % 22
y(:,4)     = sumnan(detail07(:,25:36)')';       % 23
y(:,5)     = sumnan(detail07(:,37:274)')';      % 31G
y(:,6)     = detail07(:,275);                   % 42
y(:,7)     = sumnan(detail07(:,276:279)')';     % 44RT
y(:,8)     = sumnan(detail07(:,280:288)')';     % 48TW
y(:,9)     = sumnan(detail07(:,289:303)')';     % 51
y(:,10)    = sumnan(detail07(:,304:316)')';     % FIRE
y(:,11)    = sumnan(detail07(:,317:340)')';     % PROF
y(:,12)    = sumnan(detail07(:,341:356)')';     % 6
y(:,13)    = sumnan(detail07(:,357:368)')';     % 7
y(:,14)    = sumnan(detail07(:,369:380)')';     % 81
y(:,15)    = sumnan(detail07(:,381:389)')';     % G
y          = y(:,1:15);

%
% Store intermediate purchases for certain commodities
%    * Computer systems design services
%    * Architectural, engineering, and related services
%    * Specialized design services
%    * Management consulting services
%    * Environmental and other technical consulting services
%    * Advertising, public relations, and related services
%    * Marketing research and all other misc. prof., sci., tech. services


comdes     = y(318,:);  
arch       = y(321,:);
spdes      = y(322,:);
mancon     = y(323,:);
envcon     = y(324,:);
ads        = y(326,:);
mktres     = y(327,:);
intang     = comdes+arch+spdes+mancon+envcon+ads+mktres;


%
% Combine commodities to construct finer investment splits
%

z          = zeros(15,20);
z(1,:)     = sumnan(detail07(1:13,391:410));        % 11
z(2,:)     = sumnan(detail07(14:21,391:410));       % 21
z(3,:)     = sumnan(detail07(22:24,391:410));       % 22
z(4,:)     = sumnan(detail07(25:36,391:410));       % 23
z(5,:)     = sumnan(detail07(37:273,391:410));      % 31G
z(6,:)     = detail07(274,391:410);                 % 42
z(7,:)     = sumnan(detail07(275:278,391:410));     % 44RT
z(8,:)     = sumnan(detail07(279:287,391:410));     % 48TW
z(9,:)     = sumnan(detail07(288:302,391:410));     % 51
z(10,:)    = sumnan(detail07(303:315,391:410));     % FIRE
z(11,:)    = sumnan(detail07(316:339,391:410));     % PROF
z(12,:)    = sumnan(detail07(340:355,391:410));     % 6
z(13,:)    = sumnan(detail07(356:367,391:410));     % 7
z(14,:)    = sumnan(detail07(368:379,391:410));     % 81
z(15,:)    = sumnan(detail07(380:385,391:410));     % G

Xfs        = z(:,2);
Xfe        = z(:,3);
Xfi        = z(:,4);
Xfr        = z(:,5);
Xf2        = sumnan([Xfs,Xfe,Xfi,Xfr]')';
Gdefc      = z(:,9);
Gdefs      = z(:,10);
Gdefe      = z(:,11);
Gdefi      = z(:,12);
Gndefc     = z(:,13);
Gndefs     = z(:,14);
Gndefe     = z(:,15);
Gndefi     = z(:,16);
Gslc       = z(:,17);
Gsls       = z(:,18);
Gsle       = z(:,19);
Gsli       = z(:,20);
Gc         = sumnan([Gdefc,Gndefc,Gslc]')';
Gs         = sumnan([Gdefs,Gndefs,Gsls]')';
Ge         = sumnan([Gdefe,Gndefe,Gsle]')';
Gi         = sumnan([Gdefi,Gndefi,Gsli]')';
Gx         = sumnan([Gs,Ge,Gi]')';  %G=Gc+Gx

%
% Construct capital flow tables (CFT) for nonresidential investment
% (using detailed investment data available for 63 industries, 113 years)
%

load detailinv.dat
xx         = zeros(100,63);
xx(:)      = detailinv(:,107);

%
% Group investments into 3 types and industries into 14 major categories
%
xe07       = zeros(39,14);
xs07       = zeros(32,14);
xi07       = zeros(25,14);
for i=1:63;
  if (i>=1 & i<=2); j=1; end;
  if (i>=3 & i<=5); j=2; end;
  if (i==6); j=3; end;
  if (i==7); j=4; end;
  if (i>=8 & i<=26); j=5; end;
  if (i==27); j=6; end;
  if (i==28); j=7; end;
  if (i>=29 & i<=36); j=8; end;
  if (i>=37 & i<=40); j=9; end;
  if (i>=41 & i<=47); j=10; end;
  if (i>=48 & i<=53); j=11; end;
  if (i>=54 & i<=58); j=12; end;
  if (i>=59 & i<=62); j=13; end;
  if (i==63); j=14; end;
  xe07(:,j)= xe07(:,j)+xx(3:41,i);
  xs07(:,j)= xs07(:,j)+xx(43:74,i);
  xi07(:,j)= xi07(:,j)+xx(76:100,i);
end;


%
% Initialize the capital flow tables
%
cfte            = zeros(15);
cfts            = zeros(15);
cftibea         = zeros(15);
cfti            = zeros(15);

%
% Structures:
%   * Mining structures assigned to commodity mining (21) 
%   * Everything else to construction (23)
%
cfts(2,1:14)    = sum(xs07(16:17,:));
cfts(4,1:14)    = sum(xs07([1:15,18:32],:));
cfts(2,:)       = Xfs(2)/sumnan(cfts(2,:)') *cfts(2,:);
cfts(4,:)       = Xfs(4)/sumnan(cfts(4,:)') *cfts(4,:);
cfts(5,5)       = Xfs(5);
cfts(6,6)       = Xfs(6);
cfts(8,8)       = Xfs(8);
cfts(10,10)     = Xfs(10);
cfts(:,10)      = sumnan([cfts(:,10),Xfr]')';
cfts(:,15)      = sumnan([cfts(:,15),Gs]')';

%
% Check totals against total industry FA investments in structures
%

load fa_ind_s.dat
load fa_gov.dat
fainvs          = [fa_ind_s([3,6,10:12,34:36,45],61);
                   sum(fa_ind_s([50,56],61));
                   sum(fa_ind_s([59,63,64],61));
                   sum(fa_ind_s([67,68],61));
                   sum(fa_ind_s([73,76],61));
                   fa_ind_s(79,61)
                   fa_gov(4,11)]*1000;
%
% IPP (BEA concept):
%   * Prepackaged software and artistic originals assigned to information (51)
%   % Custom and own software and all RDs to PROF 
%
cftibea(9,1:14) = xi07(1,:)+sum(xi07(21:25,:));
cftibea(11,1:14)= sum(xi07(2:20,:));
cftibea(:,15)   = sumnan([cftibea(:,15),Gi]')';

%
% IPP (more theoretically-based concept):
%   * Prepackaged software and artistic originals assigned to information (51)
%   % Custom and own software and all RDs to investing industry
%   % Defense to 31G and G, nondefense in same proportions as private, S&L to 6
%

defrd           = fa_gov(38,11)*1000;
ndefrd          = fa_gov(55,11)*1000;
slrd            = fa_gov(75,11)*1000;
cfti(9,1:14)    = xi07(1,:)+sum(xi07(21:25,:));
cfti(:,15)      = sumnan([cfti(:,15),Gi]')';
cfti(11,15)     = sumnan([cfti(11,15),-defrd,-ndefrd,-slrd]');
tem             = sum(xi07(2:20,1:14));
for i=1:14;
  cfti(i,i)     = sumnan([cfti(i,i),tem(i)]');
end;
cfti(5,5)       = cfti(5,5)+.5*defrd;
cfti(15,15)     = cfti(15,15)+.5*defrd;
cfti(12,12)     = cfti(12,12)+slrd;
tem             = ndefrd*diag(cfti)/sumnan(diag(cfti));
for i=1:15;
  cfti(i,i)     = sumnan([cfti(i,i),tem(i)]');
end;
%
% Add to cfti diagonal for FIRE based on OES wages
%
cfti(10,10)     = cfti(10,10)+255318.45;

%
% Check totals against total industry FA investments in IP
%

load fa_ind_i.dat
load fa_gov.dat
fainvi          = [fa_ind_i([3,6,10:12,34:36,45],61);
                   sum(fa_ind_i([50,56],61));
                   sum(fa_ind_i([59,63,64],61));
                   sum(fa_ind_i([67,68],61));
                   sum(fa_ind_i([73,76],61));
                   fa_ind_i(79,61)
                   fa_gov(19,11)]*1000;

%
% Equipment:
%   * First, assign commodity codes from PEQ bridge file
%   * Then distribute equipment investments to industries
%
load peqbridge07.dat
vv         = peqbridge07(1:94,:);
ccode      = vv(:,1);
itype      = -99*ones(94,1);
for i=1:94;
  if (ccode(i)>-998);
    if ((ccode(i)>=210  & ccode(i)<=219) |  ...
        (ccode(i)>=2100 & ccode(i)<=2199)) itype(i)=2; end;
    if ((ccode(i)>=300  & ccode(i)<=399) |  ...
        (ccode(i)>=3000 & ccode(i)<=3999)) itype(i)=3; end;
    if ((ccode(i)>=510  & ccode(i)<=519) |  ...
        (ccode(i)>=5100 & ccode(i)<=5199)) itype(i)=4; end;
    if ((ccode(i)>=540  & ccode(i)<=549) |  ...
        (ccode(i)>=5400 & ccode(i)<=5499)) itype(i)=5; end;
  else
    itype(i) = 0;
  end;
end;
if any(itype<0);
  disp('Warning: check mapping of commodity codes for PEQ bridge')
end;
eqsum           = sum(xe07(1:39,:));
vvsum           = sum(vv);
vvtot           = vvsum(6);
cfte(2,1:14)    = sum(vv(itype==2,2))/vvtot*eqsum;
cfte(5,1:14)    = sum(vv((itype==3 | itype==0),2))/vvtot*eqsum;
cfte(6,1:14)    = vvsum(4)/vvtot *eqsum;
cfte(7,1:14)    = vvsum(5)/vvtot *eqsum;
cfte(8,1:14)    = vvsum(3)/vvtot *eqsum;
cfte(9,1:14)    = sum(vv(itype==4,2))/vvtot*eqsum;
cfte(11,1:14)   = sum(vv(itype==5,2))/vvtot*eqsum;
cfte(:,15)      = sumnan([cfte(:,15),Ge]')';

%
% Check totals against total industry FA investments in equipment
%

load fa_ind_e.dat
load fa_gov.dat
fainve          = [fa_ind_e([3,6,10:12,34:36,45],61);
                   sum(fa_ind_e([50,56],61));
                   sum(fa_ind_e([59,63,64],61));
                   sum(fa_ind_e([67,68],61));
                   sum(fa_ind_e([73,76],61));
                   fa_ind_e(79,61)
                   fa_gov(3,11)]*1000;


%
% Make further adjustments in light of theory
%   * Adjust matrix of intermediate purchases (reducing PROF)
%

Int             = Intbea;
RK              = RKbea;
Int(11,:)       = sumnan([Int(11,:);-intang(1:15)])';
cfti(11,:)      = sumnan([cfti(11,:);intang(1:15)])';
RK              = sumnan([RK;intang(1:15)]);

%
%   * Add consumer durables and inventories to CFTs
%

load detailcd.dat
load detailcdd.dat
load detailcdk.dat
cdeq            = sum(detailcd([3:15,17:23,25:27,29:30],83));
cdsoft          = detailcd(16,83);
cdbook          = sum(detailcd([24,28],83));
cdeqdep         = sum(detailcdd([3:15,17:23,25:27,29:30],83));
cdsoftdep       = detailcdd(16,83);
cdbookdep       = sum(detailcdd([24,28],83));
cdeqk           = sum(detailcdk([3:15,17:23,25:27,29:30],83));
cdsoftk         = detailcdk(16,83);
cdbookk         = sum(detailcdk([24,28],83));

cfte(5,10)      = cfte(5,10)+cdeq;
cfti(9,10)      = cfti(9,10)+cdsoft+cdbook;
C               = sumnan([Cbea,Gc]')';
C(5)            = C(5)-cdeq+cdeqdep+.04*cdeqk;
RK(5)           = RK(5)+cdeqdep+.04*cdeqk;
RK(9)           = RK(9)+cdsoftdep+cdbookdep+.04*(cdsoftk+cdbookk);
for i=1:15;    
 cfte(i,i)      = sumnan([cfte(i,i);Xinv(i)]);
 if isnan(C(i));
   C(i)         = 0;
 end;
end;

%
%   * Reassign spending on used goods and noncomparable imports
%

Int(14,:)       = Int(14,:)+Iused+Inonc;
C(14)           = sumnan([C(14),Fused,Fnonc]');

%
%   * Distribute net exports to domestic categories
%
NX              = sumnan([EX,-IM]')';
for i=1:15;
  if (i==2);
    Int(2,:)    = Int(2,:)+NX(2)*Int(2,:)/sumnan(Int(2,:)');
    RK          = RK-NX(2)*Int(2,:)/sumnan(Int(2,:)');
  else;
    if ~isnan(NX(i))
      tem       = sumnan(Int(i,:)'+cfte(i,:)')+C(i);
      Int(i,:)  = Int(i,:)+NX(i)*Int(i,:)/tem;
      cfte(i,:) = cfte(i,:)+NX(i)*cfte(i,:)/tem;
      C(i)      = C(i)+NX(i)*C(i)/tem;
      RK        = RK-NX(i)*Int(i,:)/tem;
    end;
  end;
end;
NX              = zeros(15,1);

%
% Add everything up
%
cftt            = cfte+cfts;
Xtang           = sum(cftt')';
Xintang         = sum(cfti')';
for i=1:15;
  if isnan(C(i)); 
    C(i)        = 0;
  end;
  if isnan(Xtang(i)); 
    Xtang(i)      = 0;
  end;
  if isnan(Xintang(i)); 
    Xintang(i)      = 0;
  end;
  if isnan(NX(i)); 
    NX(i)       = 0;
  end;
end;
X               = Xtang+Xintang;
IntO            = sum(Int')';
FU              = C+X+NX;
CO              = FU+IntO;
IntI            = sum(Int);
VA              = WL+TX+RK;
IO              = VA+IntI;
GDP             = sum(VA);
GO              = sum(IO);

capshr          = RK./(IO-TX);
labshr          = WL./(IO-TX);
gamma           = Int./(ones(15,1)*(IO-TX));
zeta            = cftt./(ones(15,1)*sum(cftt));
nu              = cfti./(ones(15,1)*sum(cfti));
omega           = max(C/sum(C),.0001);
omega           = omega/sum(omega);
cftratio        = sum(cfti)'./sum(cfte+cfts)';

save gamma.dat gamma -ascii
save labshr.dat labshr -ascii
save capshr.dat capshr -ascii
save zeta.dat zeta -ascii
save nu.dat nu -ascii
save omega.dat omega -ascii
save cftratio.dat cftratio -ascii
