%MLE_Inputs    Loads in U.S. data needed for the MLE estimation. 
%              (Current results for gross output with filled-in
%               quarterly estimates based on HP filtered data are
%               stored in fillqtr1.out and corrected series---those
%               with the HP filtering done after the original trend is
%               added back---is stored in fillqtr2.out.)
%
%              The nomenclature for industry data used here is as follows.
%              The first letter is the frequency: q for quarterly and a for
%              annual. The next two letters are the variable:
%          
%                  ni     = national income 
%                  br     = business receipts
%                  em     = employment
%
%              The next two letters are the major industry (with description
%              under NAICS/SIC)
%
%                  ag     = Agriculture, forestry, fishing, and hunting/
%                           Agriculture, forestry, and fishing
%                  mi     = Mining/
%                           Mining
%                  ut     = Utilities/
%                           Electric, gas, and sanitary services
%                  co     = Construction/
%                           Construction
%                  mf     = Manufacturing/
%                           Manufacturing
%                  wt     = Wholesale trade/
%                           Wholesale trade
%                  rt     = Retail trade/
%                           Retail trade
%                  tw     = Transportation and warehousing/
%                           Transportation
%                  in     = Information/
%                           Communication
%                  fr     = Finance, insurance, real estate, rental, and
%                           leasing/Finance, insurance, and real estate
%                  pb     = Professional and business services/
%                           Share of Services
%                  eh     = Educational services, health care, and social 
%                           assistance/Share of Services
%                  ae     = Arts, entertainment, recreation, accomodation, 
%                           and food services/Share of Services
%                  os     = Other serivces/
%                           Share of Services
%                  gv     = Government/
%                           Government
%
%              If it is a 3-digit minor industry, the next set of letters are:
%
%                  fm     = Fabricated metal
%                  ma     = Machinery
%                  ce     = Computer and electronic products
%                  ee     = Electrical equipment and appliances
%                  te     = Transportation equipment
%                  ch     = Chemicals and allied products
%                  pc     = Petroleum and coal
%                  om     = Other manufacturing
%                  pu     = Publishing
%                  mp     = Motion picture and sound recording
%                  bt     = Brodcasting and telecommunications
%                  fi     = Finance and insurance
%                  pt     = Professional and technical services
%                  mc     = Management of companies
%                  os     = Other service
%                  oi     = Other industries
%
%              If it is 4-digit minor industry, the final set of letters are:
%
%                  oa     = Ordnance and accesories
%                  am     = Agricultural, mining and construction machinery
%                  im     = Industrial machinery
%                  cm     = Commercial and service machinery
%                  et     = Engine, turbine, and power transmission machinery
%                  om     = Other general purpose machinery
%                  cp     = Computer and peripheral equipment
%                  ce     = Communications equipment
%                  se     = Semiconductors and electronic components
%                  ei     = Electronic instruments
%                  mm     = Magnetic media and reproduction
%                  ee     = Electrical equipment
%                  oe     = Other electrical equipment
%                  mv     = Motor vehicles
%                  mb     = Motor vehicle bodies and trailers
%                  mp     = Motor vehicle parts
%                  ap     = Aerospace products and parts
%                  sb     = Ship and boat building
%                  ot     = Other transportation equipment
%                  me     = Medical equipment and supplies
%                  pr     = Petroleum refineries
%                  bc     = Basic chemicals
%                  rr     = Resin, rubber, and artifical fibers
%                  ac     = Agricultural chemicals
%                  pm     = Pharmaceuticals and medicines
%                  pc     = Paints, coatings and adhesives
%                  sc     = Soaps, cleaning compounds, and toiletries
%                  oc     = Other chemical products and preparations
%                  nb     = Newspapers, books, directories
%                  sf     = Software
%                  rt     = Radio and TV broadcasting
%                  cp     = Cable and other subscription programming
%                  wi     = Wired telecommunications carriers
%                  wl     = Wireless telecommunications carriers
%                  ob     = Other telecommunications
%                  dp     = Data processing, hosting, and related services
%                  oi     = Other information services
%                  dc     = Depository credit intermediation
%                  nc     = Nondepository credit intermediation
%                  ci     = Activities related to credit intermediation
%                  cc     = Securities and commodity contracts brokerage etc
%                  of     = Other financial investment activities
%                  ic     = Insurance carriers
%                  ls     = Legal services
%                  ab     = Accounting and bookkeeping
%                  ea     = Engineering and architectural
%                  cs     = Computer systems design
%                  mt     = Management and technical consulting             
%                  rd     = Scientific R&D services
%                  ad     = Advertising
%                  op     = Other professional and technical services
%                  fs     = Facilities support
%                  es     = Education support
%                  cr     = Computer and office machine repair
%                  og     = Oil and gas extraction
%                  pg     = Electric power generation
%                  wt     = Wholesale trade, professional and commercial equip.
%                  tw     = Transportation and warehousing, pipeline trans.
%
%
%              See README files in the relevant subdirectories for details 
%              on the original sources and constructed series.

%              Ellen McGrattan, 8-27-15
%              Revised, ERM, 8-30-15

clear

%======================
% BEA Population data
%======================
%
%   * Population, midperiod (thousands)

load ./rawdata/NIPA/NIPAhist/table2_1Q.dat
load ./rawdata/NIPA/NIPAhist/table2_1A.dat
popq          = table2_1Q(42,1:88);
popa          = table2_1A(41,19:40);
load ./rawdata/NIPA/NIPA/table2_1Q.dat
load ./rawdata/NIPA/NIPA/table2_1A.dat
popq          = [popq,table2_1Q(42,1:184)]';
popa          = [popa,table2_1A(41,1:46)]';
tpopq         = ptime(1947.1,272,4);
tpopa         = [1947:2014]';


%===========================
% BEA National Account data
%===========================
%
%   * GDP deflator (2009=100)
%   * Nominal GDP ($billion)
%   * Real GDP ($billion)
%
%   All series: 1947-2014, quarterly and annual
%

load ./rawdata/NIPA/NIPAhist/table1_1_9Q.dat
load ./rawdata/NIPA/NIPAhist/table1_1_5Q.dat
load ./rawdata/NIPA/NIPAhist/table1_1_6Q.dat
defq          = table1_1_9Q(3,1:88);
ngdpq         = table1_1_5Q(3,1:88);
rgdpq         = table1_1_6Q(3,1:88);

load ./rawdata/NIPA/NIPA/table1_1_9Q.dat
load ./rawdata/NIPA/NIPA/table1_1_5Q.dat
load ./rawdata/NIPA/NIPA/table1_1_6Q.dat
defq          = [defq,table1_1_9Q(3,1:184)]';
ngdpq         = [ngdpq,table1_1_5Q(3,1:184)]';
rgdpq         = [rgdpq,table1_1_6Q(3,1:184)]';

load ./rawdata/NIPA/NIPAhist/table1_1_9A.dat
load ./rawdata/NIPA/NIPAhist/table1_1_5A.dat
load ./rawdata/NIPA/NIPAhist/table1_1_6A.dat
defa          = table1_1_9A(2,19:40);
ngdpa         = table1_1_5A(2,19:40);
rgdpa         = table1_1_6A(2,19:40);

load ./rawdata/NIPA/NIPA/table1_1_9A.dat
load ./rawdata/NIPA/NIPA/table1_1_5A.dat
load ./rawdata/NIPA/NIPA/table1_1_6A.dat
defa          = [defa,table1_1_9A(2,1:46)]';
ngdpa         = [ngdpa,table1_1_5A(2,1:46)]';
rgdpa         = [rgdpa,table1_1_6A(2,1:46)]';
tgdpq         = ptime(1947.1,272,4);
tgdpa         = [1947:2014]';


%===========================
% BEA National Incomes data
%===========================
%
%   * National incomes (wo/ ccadj) by NAICs/SIC industry ($billion)
%
%   All series: 1948-2014, quarterly and annual (NaNs added for 1947)
%
load ./rawdata/NIPA/NIPA/table6_1DQ.dat
qni           = table6_1DQ(3,1:56);
qniag         = table6_1DQ(6,1:56);
qnimi         = table6_1DQ(7,1:56);
qniut         = table6_1DQ(8,1:56);
qnico         = table6_1DQ(9,1:56);
qnimf         = table6_1DQ(10,1:56);
qniwt         = table6_1DQ(13,1:56);
qnirt         = table6_1DQ(14,1:56);
qnitw         = table6_1DQ(15,1:56);
qniin         = table6_1DQ(16,1:56);
qnifr         = table6_1DQ(17,1:56);
qnipb         = table6_1DQ(18,1:56);
qnieh         = table6_1DQ(19,1:56);
qniae         = table6_1DQ(20,1:56);
qnios         = table6_1DQ(21,1:56);
qnigv         = table6_1DQ(22,1:56);

load ./rawdata/NIPA/NIPA/table6_1CQ.dat
load ./rawdata/NIPA/NIPA/table6_1BQ.dat
tab6          = [table6_1BQ(:,1:72),table6_1CQ(:,1:56)];

load ./rawdata/NIPA/NIPAhist/table6_1BQ.dat
tab6          = [table6_1BQ(:,1:84),tab6];
grate         = qni(1)/tab6(3,212);
qni           = [tab6(3,:),qni]';
qniag         = [tab6(6,:) /tab6(6,212) *qniag(1)/grate,qniag]';
qnimi         = [tab6(7,:) /tab6(7,212) *qnimi(1)/grate,qnimi]';
qniut         = [tab6(15,:)/tab6(15,212)*qniut(1)/grate,qniut]';
qnico         = [tab6(8,:) /tab6(8,212) *qnico(1)/grate,qnico]';
qnimf         = [tab6(9,:) /tab6(9,212) *qnimf(1)/grate,qnimf]';
qniwt         = [tab6(16,:)/tab6(16,212)*qniwt(1)/grate,qniwt]';
qnirt         = [tab6(17,:)/tab6(17,212)*qnirt(1)/grate,qnirt]';
qnitw         = [tab6(13,:)/tab6(13,212)*qnitw(1)/grate,qnitw]';
qniin         = [tab6(14,:)/tab6(14,212)*qniin(1)/grate,qniin]';
qnifr         = [tab6(18,:)/tab6(18,212)*qnifr(1)/grate,qnifr]';
qnipb         = [tab6(19,:)/tab6(19,212)*qnipb(1)/grate,qnipb]';
qnieh         = [tab6(19,:)/tab6(19,212)*qnieh(1)/grate,qnieh]';
qniae         = [tab6(19,:)/tab6(19,212)*qniae(1)/grate,qniae]';
qnios         = [tab6(19,:)/tab6(19,212)*qnios(1)/grate,qnios]';
qnigv         = [tab6(20,:)/tab6(20,212)*qnigv(1)/grate,qnigv]';
qni15         = [qniag,qnimi,qniut,qnico,qnimf,qniwt,qnirt,qnitw, ...
                 qniin,qnifr,qnipb,qnieh,qniae,qnios,qnigv];

ani           = zeros(4,67); ani(:)   = qni;   ani   = mean(ani)';
aniag         = zeros(4,67); aniag(:) = qniag; aniag = mean(aniag)';
animi         = zeros(4,67); animi(:) = qnimi; animi = mean(animi)';
aniut         = zeros(4,67); aniut(:) = qniut; aniut = mean(aniut)';
anico         = zeros(4,67); anico(:) = qnico; anico = mean(anico)';
animf         = zeros(4,67); animf(:) = qnimf; animf = mean(animf)';
aniwt         = zeros(4,67); aniwt(:) = qniwt; aniwt = mean(aniwt)';
anirt         = zeros(4,67); anirt(:) = qnirt; anirt = mean(anirt)';
anitw         = zeros(4,67); anitw(:) = qnitw; anitw = mean(anitw)';
aniin         = zeros(4,67); aniin(:) = qniin; aniin = mean(aniin)';
anifr         = zeros(4,67); anifr(:) = qnifr; anifr = mean(anifr)';
anipb         = zeros(4,67); anipb(:) = qnipb; anipb = mean(anipb)';
anieh         = zeros(4,67); anieh(:) = qnieh; anieh = mean(anieh)';
aniae         = zeros(4,67); aniae(:) = qniae; aniae = mean(aniae)';
anios         = zeros(4,67); anios(:) = qnios; anios = mean(anios)';
anigv         = zeros(4,67); anigv(:) = qnigv; anigv = mean(anigv)';
ani15         = [aniag,animi,aniut,anico,animf,aniwt,anirt,anitw, ...
                 aniin,anifr,anipb,anieh,aniae,anios,anigv];
ani2          = sum(ani15')';
tniq          = ptime(1948.1,268,4);
tnia          = [1948:2014]';


%=========================================
% SIC/NAICS income ratio for 1997 to 1998
%=========================================
% 
%   * Based on BEA national incomes
%   
%
t97           = find(tnia==1997);
t98           = find(tnia==1998);

rat           = ani(t97)/ani(t98);
rat15         = ani15(t97,:)./ani15(t98,:);

%====================
% BLS Population data
%====================
%
%   * CPS noninstitutional population, ages 16-64
%

load ./rawdata/Employment/CPU/pop1664.dat

p16q          = [NaN;NaN;pop1664(1:270,2)];
tp16q         = [1947:2014]';


%====================
% BLS Employment data
%====================
%
%   * CES and EES, major and some minor industries (thousands)
%  
%   All series: 1947-2014, quarterly and annual


load ./rawdata/Employment/BLS/ce_maj.dat
load ./rawdata/Employment/BLS/ee_maj.dat
load ./rawdata/Employment/Ramey/employment.dat

qem           = zeros(272,1);
aem           = zeros(68,1);
qem15         = zeros(272,15);
aem15         = zeros(68,15);
qem15(:,1)    = employment(1:272,13)*1000;
tem           = zeros(4,68);
tem(:)        = qem15(:,1);
aem15(:,1)    = mean(tem)';

for i=2:15;
  ce          = zeros(3,272);
  ce(:)       = ce_maj(1:816,i-1);
  mce         = mean(ce)';
  if (i==3);
    mce(1:205)= mean(ee_maj(1:205,:)')'*mce(205)/mean(ee_maj(205,:));
    mce(1:37) = ee_maj(1:37,1)*mce(37)/ee_maj(37,1);
  end;
  qem15(:,i)  = mce;
  tem(:)      = qem15(:,i);
  aem15(:,i)  = mean(tem)';
end;
ce(:)         = ce_maj(1:816,15);
qem           = mean(ce)'+qem15(:,1);
tem(:)        = qem;
aem           = mean(tem)';
qemag         = qem15(:,1);
qemmi         = qem15(:,2);
qemut         = qem15(:,3);
qemco         = qem15(:,4);
qemmf         = qem15(:,5);
qemwt         = qem15(:,6);
qemrt         = qem15(:,7);
qemtw         = qem15(:,8);
qemin         = qem15(:,9);
qemfr         = qem15(:,10);
qempb         = qem15(:,11);
qemeh         = qem15(:,12);
qemae         = qem15(:,13);
qemos         = qem15(:,14);
qemgv         = qem15(:,15);
aemag         = aem15(:,1);
aemmi         = aem15(:,2);
aemut         = aem15(:,3);
aemco         = aem15(:,4);
aemmf         = aem15(:,5);
aemwt         = aem15(:,6);
aemrt         = aem15(:,7);
aemtw         = aem15(:,8);
aemin         = aem15(:,9);
aemfr         = aem15(:,10);
aempb         = aem15(:,11);
aemae         = aem15(:,12);
aemeh         = aem15(:,13);
aemos         = aem15(:,14);
aemgv         = aem15(:,15);


load ./rawdata/Employment/BLS/ce_min.dat
load ./rawdata/Employment/BLS/ee_min.dat
qem60         = zeros(272,15);
aem60         = zeros(68,15);
for i=1:60;
  ce          = zeros(3,272);
  ee          = zeros(3,272);
  ce(:)       = ce_min(1:816,i);
  ee(:)       = ee_min(1:816,i);
  mce         = mean(ce)';
  mee         = mean(ee)';
  j           = find(~isnan(mce));
  j           = j(1);
  mce(1:j)    = mee(1:j)*mce(j)/mee(j);
  qem60(:,i)  = mce;
  tem(:)      = mce;
  aem60(:,i)  = mean(tem)';
end;

qemmffmoa     = qem60(:,1);
qemmfmaam     = qem60(:,2);
qemmfmaim     = qem60(:,3);
qemmfmacm     = qem60(:,4);
qemmfmaet     = qem60(:,5);
qemmfmaom     = qem60(:,6);
qemmfcecp     = qem60(:,7);
qemmfcece     = qem60(:,8);
qemmfcese     = qem60(:,9);
qemmfceei     = qem60(:,10);
qemmfcemm     = qem60(:,11);
qemmfeeee     = qem60(:,12);
qemmfeeoe     = qem60(:,13);
qemmftemv     = qem60(:,14);
qemmftemb     = qem60(:,15);
qemmftemp     = qem60(:,16);
qemmfteap     = qem60(:,17);
qemmftesb     = qem60(:,18);
qemmfteot     = qem60(:,19);
qemmfomme     = qem60(:,20);
qemmfpcpr     = qem60(:,21);
qemmfchbc     = qem60(:,22);
qemmfchrr     = qem60(:,23);
qemmfchac     = qem60(:,24);
qemmfchpm     = qem60(:,25);
qemmfchpc     = qem60(:,26);
qemmfchsc     = qem60(:,27);
qemmfchoc     = qem60(:,28);
qeminpunb     = qem60(:,29);
qeminpusf     = qem60(:,30);
qeminmp       = qem60(:,31);
qeminbtrt     = qem60(:,32);
qeminbtcp     = qem60(:,33);
qeminbtwi     = qem60(:,34);
qeminbtwl     = qem60(:,35);
qeminbtob     = qem60(:,36);
qemindp       = qem60(:,37);
qeminoi       = qem60(:,38);
qemfrfidc     = qem60(:,39);
qemfrfinc     = qem60(:,40);
qemfrfici     = qem60(:,41);
qemfrficc     = qem60(:,42);
qemfrfiof     = qem60(:,43);
qemfrfiic     = qem60(:,44);
qempbptls     = qem60(:,45);
qempbptab     = qem60(:,46);
qempbptea     = qem60(:,47);
qempbptcs     = qem60(:,48);
qempbptmt     = qem60(:,49);
qempbptrd     = qem60(:,50);
qempbptad     = qem60(:,51);
qempbptop     = qem60(:,52);
qempbmc       = qem60(:,53);
qempbadfs     = qem60(:,54);
qemehedes     = qem60(:,55);
qemoscr       = qem60(:,56);
qemmiog       = qem60(:,57);
qemutpg       = qem60(:,58);
qemwtpe       = qem60(:,59);
qemtwpt       = qem60(:,60);

qemmfch       = zeros(272,1);
ii            = 101:272;
qemmfch(ii)   = sum(qem60(ii,22:27)')';
ii            = 45:100;
qemmfch(ii)   = sum(qem60(ii,[23,25:27])')';
qemmfch(ii)   = qemmfch(ii)*qemmfch(101)/qemmfch(100);
ii            = 1:44;
qemmfch(ii)   = qem60(ii,25);
qemmfch(ii)   = qemmfch(ii)*qemmfch(45)/qemmfch(44);

qeminbt       = zeros(272,1);
ii            = 165:272;
qeminbt(ii)   = sum(qem60(ii,32:37)')';
ii            = 101:164;
qeminbt(ii)   = sum(qem60(ii,[32,34:37])')';
qeminbt(ii)   = qeminbt(ii)*qeminbt(165)/qeminbt(164);
ii            = 45:100;
qeminbt(ii)   = sum(qem60(ii,[32,34:36])')';
qeminbt(ii)   = qeminbt(ii)*qeminbt(101)/qeminbt(100);
ii            = 1:44;
qeminbt(ii)   = sum(qem60(ii,[34:36])')';
qeminbt(ii)   = qeminbt(ii)*qeminbt(45)/qeminbt(44);

qeminpusf     = zeros(272,1);
ii            = 101:272;
qeminpusf(ii) = qem60(ii,30);
ii            = 1:100;
qeminpusf(ii) = qemin(ii)*qeminpusf(101)/qemin(101);


qemmfce       = zeros(272,1);
ii            = 165:272;
qemmfce(ii)   = sum(qem60(ii,7:11)')';
ii            = 45:164;
qemmfce(ii)   = sum(qem60(ii,[7:9,11])')';
qemmfce(ii)   = qemmfce(ii)*qemmfce(165)/qemmfce(164);
ii            = 1:44;
qemmfce(ii)   = qemmf(ii)*qemmfce(45)/qemmf(45);

qemfrfi       = zeros(272,1);
ii            = 165:272;
qemfrfi(ii)   = sum(qem60(ii,39:44)')';
ii            = 141:164;
qemfrfi(ii)   = sum(qem60(ii,[39,41,42,44])')';
qemfrfi(ii)   = qemfrfi(ii)*qemfrfi(165)/qemfrfi(164);
ii            = 101:140;
qemfrfi(ii)   = sum(qem60(ii,[39,42,44])')';
qemfrfi(ii)   = qemfrfi(ii)*qemfrfi(141)/qemfrfi(140);
ii            = 45:100;
qemfrfi(ii)   = qem60(ii,44);
qemfrfi(ii)   = qemfrfi(ii)*qemfrfi(101)/qemfrfi(100);
ii            = 1:44;
qemfrfi(ii)   = qemfr(ii)*qemfrfi(45)/qemfr(45);


aemmffmoa     = aem60(:,1);
aemmfmaam     = aem60(:,2);
aemmfmaim     = aem60(:,3);
aemmfmacm     = aem60(:,4);
aemmfmaet     = aem60(:,5);
aemmfmaom     = aem60(:,6);
aemmfcecp     = aem60(:,7);
aemmfcece     = aem60(:,8);
aemmfcese     = aem60(:,9);
aemmfceei     = aem60(:,10);
aemmfcemm     = aem60(:,11);
aemmfeeee     = aem60(:,12);
aemmfeeoe     = aem60(:,13);
aemmftemv     = aem60(:,14);
aemmftemb     = aem60(:,15);
aemmftemp     = aem60(:,16);
aemmfteap     = aem60(:,17);
aemmftesb     = aem60(:,18);
aemmfteot     = aem60(:,19);
aemmfomme     = aem60(:,20);
aemmfpcpr     = aem60(:,21);
aemmfchbc     = aem60(:,22);
aemmfchrr     = aem60(:,23);
aemmfchac     = aem60(:,24);
aemmfchpm     = aem60(:,25);
aemmfchpc     = aem60(:,26);
aemmfchsc     = aem60(:,27);
aemmfchoc     = aem60(:,28);
aeminpunb     = aem60(:,29);
aeminpusf     = aem60(:,30);
aeminmp       = aem60(:,31);
aeminbtrt     = aem60(:,32);
aeminbtcp     = aem60(:,33);
aeminbtwi     = aem60(:,34);
aeminbtwl     = aem60(:,35);
aeminbtob     = aem60(:,36);
aemindp       = aem60(:,37);
aeminoi       = aem60(:,38);
aemfrfidc     = aem60(:,39);
aemfrfinc     = aem60(:,40);
aemfrfici     = aem60(:,41);
aemfrficc     = aem60(:,42);
aemfrfiof     = aem60(:,43);
aemfrfiic     = aem60(:,44);
aempbptls     = aem60(:,45);
aempbptab     = aem60(:,46);
aempbptea     = aem60(:,47);
aempbptcs     = aem60(:,48);
aempbptmt     = aem60(:,49);
aempbptrd     = aem60(:,50);
aempbptad     = aem60(:,51);
aempbptop     = aem60(:,52);
aempbmc       = aem60(:,53);
aempbadfs     = aem60(:,54);
aemehedes     = aem60(:,55);
aemoscr       = aem60(:,56);
aemmiog       = aem60(:,57);
aemutpg       = aem60(:,58);
aemwtpe       = aem60(:,59);
aemtwpt       = aem60(:,60);

temq          = ptime(1947.1,272,4);
tema          = [1947:2014]';

%
%   * Labor productivity analyses, numbers of hours and employees
%  
%   All series: 1987-2014, annual 
%

load ./rawdata/Employment/LPC/lpc_hours.dat
load ./rawdata/Employment/LPC/lpc_empl.dat
lpc_hours     = lpc_hours*10^6;
lpc_empl      = lpc_empl*10^3;
thra          = [1947:2014]';
thrq          = ptime(1947.1,272,4);
ii            = [1,8,18,28,34];
anhag         = [sum(lpc_hours(ii,1))*ones(1,40),sum(lpc_hours(ii,1:28))]';
aneag         = [sum(lpc_empl(ii,1))*ones(1,40),sum(lpc_empl(ii,1:28))]';
tem           = ann2qtr([thra,anhag./aneag],3,0); 
qhrag         = tem(:,2).*qemag;
anhmi         = [lpc_hours(43,1)*ones(1,40),lpc_hours(43,1:28)]';
anemi         = [lpc_empl(43,1)*ones(1,40),lpc_empl(43,1:28)]';
tem           = ann2qtr([thra,anhmi./anemi],3,0); 
qhrmi         = tem(:,2).*qemmi;
anhut         = [lpc_hours(55,1)*ones(1,40),lpc_hours(55,1:28)]';
aneut         = [lpc_empl(55,1)*ones(1,40),lpc_empl(55,1:28)]';
tem           = ann2qtr([thra,anhut./aneut],3,0); 
qhrut         = tem(:,2).*qemut;
ii            = [62,66,77];
anhco         = [sum(lpc_hours(ii,1))*ones(1,40),sum(lpc_hours(ii,1:28))]';
aneco         = [sum(lpc_empl(ii,1))*ones(1,40),sum(lpc_empl(ii,1:28))]';
tem           = ann2qtr([thra,anhco./aneco],3,0); 
qhrco         = tem(:,2).*qemco;
ii            = [82,102,107,113,116,122,131,142,149,151,155,164,170,179, ...
                 188,216,235,252,258,289,300];
anhmf         = [sum(lpc_hours(ii,1))*ones(1,40),sum(lpc_hours(ii,1:28))]';
anemf         = [sum(lpc_empl(ii,1))*ones(1,40),sum(lpc_empl(ii,1:28))]';
tem           = ann2qtr([thra,anhmf./anemf],3,0); 
qhrmf         = tem(:,2).*qemmf;
anhmfch       = [lpc_hours(155,1)*ones(1,40),lpc_hours(155,1:28)]';
anemfch       = [lpc_empl(155,1)*ones(1,40),lpc_empl(155,1:28)]';
tem           = ann2qtr([thra,anhmfch./anemfch],3,0); 
qhrmfch       = tem(:,2).*qemmfch;

anhmfce       = [lpc_hours(235,1)*ones(1,40),lpc_hours(235,1:28)]';
anemfce       = [lpc_empl(235,1)*ones(1,40),lpc_empl(235,1:28)]';
tem           = ann2qtr([thra,anhmfce./anemfce],3,0); 
qhrmfce       = tem(:,2).*qemmfce;

anhmfcese     = [lpc_hours(244,1)*ones(1,40),lpc_hours(244,1:28)]';
anemfcese     = [lpc_empl(244,1)*ones(1,40),lpc_empl(244,1:28)]';
tem           = ann2qtr([thra,anhmfcese./anemfcese],3,0); 
qhrmfcese     = tem(:,2).*qemmfcese;

anhmftemv     = [lpc_hours(259,1)*ones(1,40),lpc_hours(259,1:28)]';
anemftemv     = [lpc_empl(259,1)*ones(1,40),lpc_empl(259,1:28)]';
tem           = ann2qtr([thra,anhmftemv./anemftemv],3,0); 
qhrmftemv     = tem(:,2).*qemmftemv;

anhwt         = [lpc_hours(308,1)*ones(1,40),lpc_hours(308,1:28)]';
anewt         = [lpc_empl(308,1)*ones(1,40),lpc_empl(308,1:28)]';
tem           = ann2qtr([thra,anhwt./anewt],3,0); 
qhrwt         = tem(:,2).*qemwt;
anhrt         = [lpc_hours(333,1)*ones(1,40),lpc_hours(333,1:28)]';
anert         = [lpc_empl(333,1)*ones(1,40),lpc_empl(333,1:28)]';
tem           = ann2qtr([thra,anhrt./anert],3,0); 
qhrrt         = tem(:,2).*qemrt;
ii            = [447,452,456,461,473,488,496,506,518,522,529];
anhtw         = [sum(lpc_hours(ii,1))*ones(1,40),sum(lpc_hours(ii,1:28))]';
anetw         = [sum(lpc_empl(ii,1))*ones(1,40),sum(lpc_empl(ii,1:28))]';
tem           = ann2qtr([thra,anhtw./anetw],3,0); 
qhrtw         = tem(:,2).*qemtw;
ii            = [535,546,551,553,557,560,566,568,572];
anhin         = [sum(lpc_hours(ii,1))*ones(1,40),sum(lpc_hours(ii,1:28))]';
anein         = [sum(lpc_empl(ii,1))*ones(1,40),sum(lpc_empl(ii,1:28))]';
tem           = ann2qtr([thra,anhin./anein],3,0); 
qhrin         = tem(:,2).*qemin;
anhinpusf     = [lpc_hours(543,1)*ones(1,40),lpc_hours(543,1:28)]';
aneinpusf     = [lpc_empl(543,1)*ones(1,40),lpc_empl(543,1:28)]';
tem           = ann2qtr([thra,anhinpusf./aneinpusf],3,0); 
qhrinpusf     = tem(:,2).*qeminpusf;
ii            = [551,553,557,560,566,568,572];
anhinbt       = [sum(lpc_hours(ii,1))*ones(1,40),sum(lpc_hours(ii,1:28))]';
aneinbt       = [sum(lpc_empl(ii,1))*ones(1,40),sum(lpc_empl(ii,1:28))]';
tem           = ann2qtr([thra,anhinbt./aneinbt],3,0); 
qhrinbt       = tem(:,2).*qeminbt;

ii            = [574,578,584,590,593];
anhfrfi       = [sum(lpc_hours(ii,1))*ones(1,40),sum(lpc_hours(ii,1:28))]';
anefrfi       = [sum(lpc_empl(ii,1))*ones(1,40),sum(lpc_empl(ii,1:28))]';
tem           = ann2qtr([thra,anhfrfi./anefrfi],3,0); 
qhrfrfi       = tem(:,2).*qemfrfi;
ii            = [596,602,614];
anhfrre       = [sum(lpc_hours(ii,1))*ones(1,40),sum(lpc_hours(ii,1:28))]';
anefrre       = [sum(lpc_empl(ii,1))*ones(1,40),sum(lpc_empl(ii,1:28))]';
anhfr         = anhfrfi+anhfrre;
anefr         = anefrfi+anefrre;
tem           = ann2qtr([thra,anhfr./anefr],3,0); 
qhrfr         = tem(:,2).*qemfr;
ii            = [618,641,645,663];
anhpb         = [sum(lpc_hours(ii,1))*ones(1,40),sum(lpc_hours(ii,1:28))]';
anepb         = [sum(lpc_empl(ii,1))*ones(1,40),sum(lpc_empl(ii,1:28))]';
tem           = ann2qtr([thra,anhpb./anepb],3,0); 
qhrpb         = tem(:,2).*qempb;
anhpbptrd     = [lpc_hours(635,1)*ones(1,40),lpc_hours(635,1:28)]';
anepbptrd     = [lpc_empl(635,1)*ones(1,40),lpc_empl(635,1:28)]';
tem           = ann2qtr([thra,anhpbptrd./anepbptrd],3,0); 
qhrpbptrd     = tem(:,2).*qempbptrd;

anhpbptad     = [lpc_hours(636,1)*ones(1,40),lpc_hours(636,1:28)]';
anepbptad     = [lpc_empl(636,1)*ones(1,40),lpc_empl(636,1:28)]';
tem           = ann2qtr([thra,anhpbptad./anepbptad],3,0); 
qhrpbptad     = tem(:,2).*qempbptad;

anhpbptea     = [lpc_hours(626,1)*ones(1,40),lpc_hours(626,1:28)]';
anepbptea     = [lpc_empl(626,1)*ones(1,40),lpc_empl(626,1:28)]';
tem           = ann2qtr([thra,anhpbptea./anepbptea],3,0); 
qhrpbptea     = tem(:,2).*qempbptea;

ii            = [669,687,703,714,724];
anheh         = [sum(lpc_hours(ii,1))*ones(1,40),sum(lpc_hours(ii,1:28))]';
aneeh         = [sum(lpc_empl(ii,1))*ones(1,40),sum(lpc_empl(ii,1:28))]';
tem           = ann2qtr([thra,anheh./aneeh],3,0); 
qhreh         = tem(:,2).*qemeh;
ii            = [733,744,746,758];
anhae         = [sum(lpc_hours(ii,1))*ones(1,40),sum(lpc_hours(ii,1:28))]';
aneae         = [sum(lpc_empl(ii,1))*ones(1,40),sum(lpc_empl(ii,1:28))]';
tem           = ann2qtr([thra,anhae./aneae],3,0); 
qhrae         = tem(:,2).*qemae;
ii            = [777,787,801,813];
anhos         = [sum(lpc_hours(ii,1))*ones(1,40),sum(lpc_hours(ii,1:28))]';
aneos         = [sum(lpc_empl(ii,1))*ones(1,40),sum(lpc_empl(ii,1:28))]';
tem           = ann2qtr([thra,anhos./aneos],3,0); 
qhros         = tem(:,2).*qemos;
anhgv         = [lpc_hours(817,1)*ones(1,40),lpc_hours(817,1:28)]';
anegv         = [lpc_empl(817,1)*ones(1,40),lpc_empl(817,1:28)]';
tem           = ann2qtr([thra,anhgv./anegv],3,0); 
qhrgv         = tem(:,2).*qemgv;

anh15         = [anhag,anhmi,anhut,anhco,anhmf,anhwt,anhrt,anhtw, ...
                 anhin,anhfr,anhpb,anheh,anhae,anhos,anhgv];
ane15         = [aneag,anemi,aneut,aneco,anemf,anewt,anert,anetw, ...
                 anein,anefr,anepb,aneeh,aneae,aneos,anegv];
qhr15         = [qhrag,qhrmi,qhrut,qhrco,qhrmf,qhrwt,qhrrt,qhrtw, ...
                 qhrin,qhrfr,qhrpb,qhreh,qhrae,qhros,qhrgv] *1000;
qhrintan      = [qhrmfch,qhrinbt,qhrinpusf,qhrpbptad,qhrmfce, ...
                 qhrmftemv]*1000;
anh           = sum(anh15')';
ane           = sum(ane15')';
qhr           = sum(qhr15')';
qhpc          = qhr./p16q/5200;
qhpc21        = [qhr15,qhrintan,qhr]./(p16q*ones(1,22))/5200; 
qhpcext       = [qhrmfcese,qhrpbptea]./(p16q*ones(1,2))/5200;
for i=1:22;
  ldhpc(:,i)  = log(qhpc21(:,i)/mean(qhpc21(209:268,i)));
end;
ldhpc         = ldhpc(5:272,:);  % saved as ldhpc21.dat
for i=1:2;
  ldhpcext(:,i) = log(qhpcext(:,i)/mean(qhpcext(209:268,i)));
end;
ldhpcext      = ldhpcext(5:272,:);

load ./rawdata/Employment/CPU/avghrs.dat
qhpc2         = [NaN;NaN;avghrs(1:270,2)];
ldhpc2        = log(qhpc2/mean(qhpc2(209:268)));
ldhpc2        = ldhpc2(5:272);   % saved as ldhpctot.dat


%=======================
% IRS Business Tax data
%=======================
% 
%   * Nominal business receipts, 1929-2012, annual ($thousands)
%   
%

load ./rawdata/IRS/busrecpt/prop5798.dat
load ./rawdata/IRS/busrecpt/prop9812.dat
load ./rawdata/IRS/busrecpt/part5797.dat
load ./rawdata/IRS/busrecpt/part9812.dat
load ./rawdata/IRS/busrecpt/corp2997.dat
load ./rawdata/IRS/busrecpt/corp9812.dat

%
% Convert all figures to billions
%
prop5798      = [prop5798(:,1),prop5798(:,2:40)/10^6];
prop9812      = [prop9812(:,1),prop9812(:,2:45)/10^6];
part5797      = [part5797(:,1),part5797(:,2:40)/10^6];
part9812      = [part9812(:,1),part9812(:,2:42)/10^6];
corp2997      = [corp2997(:,1),corp2997(:,2:36)/10^6];
corp9812      = [corp9812(1,:);corp9812(2:42,:)/10^6];

%
% Use corporate data for 1947-2012
%
corp1         = corp2997(19:69,:);
corp2         = corp9812';

%
% Fill in NaNs 
%
corp1(1,12)   = corp1(2,12)*corp1(1,13)/corp1(2,13);
corp1(6,29)   = corp1(7,29)*corp1(6,28)/corp1(7,28);
corp1(6,30)   = corp1(6,28)-corp1(6,29);
corp1(6,20)   = .5*corp1(5,20)+.5*corp1(7,20);
corp1(6,31)   = .5*corp1(5,31)+.5*corp1(7,31);
corp1(6,18)   = sum(corp1(6,19:21));
i             = 1:16;
j             = 32:34;
corp1(i,j)    = corp1(i,25)*corp1(17,j)/corp1(17,25);
i             = 1:30;
j             = 32:36;
corp1(i,j)    = corp1(i,25)*corp1(31,j)/corp1(31,25);
corp1(22,33)  = .5*corp1(21,33)+.5*corp1(23,33);

%
% Scale FIRE between 1957 and 1958
%
i             = 1:11;
j             = 23;
corp1(i,j)    = corp1(i,j)*corp1(12,j)*corp1(12,j)/corp1(11,j)/corp1(13,j);

%
% Use scaled propreitor and partner returns for health services, 1987-1993
%
i0            = find(corp1(:,1)==1986);
j0            = find(corp1(:,1)==1994);
i1            = find(prop5798(:,1)==1986);
j1            = find(prop5798(:,1)==1994);
i2            = find(part5797(:,1)==1986);
j2            = find(part5797(:,1)==1994);
tem           = prop5798(i1:j1,32)+part5797(i2:j2,32);
ii            = i0+1:j0;
corp1(ii,32)  = tem(2:9)*corp1(j0,32)/tem(9);
ii            = 1:i0;
corp1(ii,32)  = corp1(ii,32)*corp1(i0,32)/corp1(i0-1,32);


sumce         = sum(corp1(:,12:13)')';
sumbt         = sum(corp2(:,20:23)')';
sumcd         = sum(corp2(:,30:31)')';
sumpb         = sum(corp2(:,[26,36:37])')';
summe         = sum(corp1(:,[32,34])')';
sumeh         = sum(corp2(:,38:39)')';
sumae         = sum(corp2(:,40:41)')';
sumha         = sum(corp1(:,[26,31])')';
sumte         = sum(corp1(:,10:11)')';

abr           = [corp1(:,2) *corp2(1,2) /corp1(51,2) *rat;     corp2(:,2)];
abrag         = [corp1(:,3) *corp2(1,3) /corp1(51,3) *rat15(1);corp2(:,3)];
abrmi         = [corp1(:,4) *corp2(1,4) /corp1(51,4) *rat15(2);corp2(:,4)];
abrut         = [corp1(:,5) *corp2(1,5) /corp1(51,5) *rat15(3);corp2(:,5)];
abrco         = [corp1(:,6) *corp2(1,6) /corp1(51,6) *rat15(4);corp2(:,6)];
abrmf         = [corp1(:,7) *corp2(1,7) /corp1(51,7) *rat15(5);corp2(:,7)];
abrmfch       = [corp1(:,8) *corp2(1,8) /corp1(51,8) *rat15(5);corp2(:,8)];
abrmfce       = [sumce      *corp2(1,9) /sumce(51)   *rat15(5);corp2(:,9)];
abrmfee       = [corp1(:,9) *corp2(1,10)/corp1(51,9) *rat15(5);corp2(:,10)];
abrmfte       = [sumte      *corp2(1,11)/sumte(51)   *rat15(5);corp2(:,11)];
abrmftemv     = [corp1(:,10)*corp2(1,12)/corp1(51,10)*rat15(5);corp2(:,12)];
abrwt         = [corp1(:,15)*corp2(1,14)/corp1(51,15)*rat15(6);corp2(:,14)];
abrrt         = [corp1(:,16)*corp2(1,15)/corp1(51,16)*rat15(7);corp2(:,15)];
abrtw         = [corp1(:,17)*corp2(1,16)/corp1(51,17)*rat15(8);corp2(:,16)];
abrinpu       = [corp1(:,19)*corp2(1,18)/corp1(51,19)*rat15(9);corp2(:,18)];
abrinmp       = [corp1(:,20)*corp2(1,19)/corp1(51,20)*rat15(9);corp2(:,19)];
abrinbt       = [corp1(:,21)*sumbt(1)   /corp1(51,21)*rat15(9);sumbt];
abrin         = abrinpu+abrinmp+abrinbt;
abrfrfi       = [corp1(:,23)*corp2(1,24)/corp1(51,23)*rat15(10);corp2(:,24)];
abrfrre       = [corp1(:,24)*corp2(1,25)/corp1(51,24)*rat15(10);corp2(:,25)];
abrfr         = abrfrfi+abrfrre;
abrpbptls     = [corp1(:,33)*corp2(1,27)/corp1(51,33)*rat15(11);corp2(:,27)];
abrpbptab     = [corp1(:,36)*corp2(1,28)/corp1(51,36)*rat15(11);corp2(:,28)];
abrpbptea     = [corp1(:,35)*corp2(1,29)/corp1(51,35)*rat15(11);corp2(:,29)];
abrpbptcs     = [corp1(:,30)*sumcd(1)   /corp1(51,30)*rat15(11);sumcd];
abrpbptmt     = [corp1(:,30)*corp2(1,32)/corp1(51,30)*rat15(11);corp2(:,32)];
abrpbptrd     = [corp1(:,30)*corp2(1,33)/corp1(51,30)*rat15(11);corp2(:,33)];
abrpbptad     = [corp1(:,29)*corp2(1,34)/corp1(51,29)*rat15(11);corp2(:,34)];
abrpbptop     = [corp1(:,30)*corp2(1,35)/corp1(51,30)*rat15(11);corp2(:,35)];
abrpbmc       = [corp1(:,30)*corp2(1,36)/corp1(51,30)*rat15(11);corp2(:,36)];
abrpb         = [corp1(:,28)*sumpb(1)   /corp1(51,28)*rat15(11);sumpb];
abreh         = [summe      *sumeh(1)   /summe(51)   *rat15(12);sumeh];
abrae         = [sumha      *sumae(1)   /sumha(51)   *rat15(13);sumae];
abros         = [corp1(:,27)*corp2(1,42)/corp1(51,27)*rat15(14);corp2(:,42)];
abrgv         = NaN*abr;

abr15         = [abrag,abrmi,abrut,abrco,abrmf,abrwt,abrrt,abrtw, ...
                 abrin,abrfr,abrpb,abreh,abrae,abros,abrgv];
abr2          = sum(abr15')';
tbra          = [1947:2012]';

%=================================
% BEA Gross Output and Value Added
%=================================
%
%   * Nominal data, annual, $millions
%   * Nominal data, quarterly, $billions
%
%   All series: 1947-2014
%


load ./rawdata/Industry/Annual/GO47to15.dat
load ./rawdata/Industry/Annual/GO403ind.dat
load ./rawdata/Industry/Annual/GO403intan.dat
load ./rawdata/Industry/Qtrly/GO05to15Q.dat
load ./rawdata/Industry/Qtrly/GO12to15Qdetail.dat

%
% Convert all figures to $billions
%

GO403ind      = [GO403ind(1,:);GO403ind(2:404,:)/1000];
GO403intan    = [GO403intan(1,:);GO403intan(2:108,:)/1000];
GO47to15      = GO47to15(2:102,:)/1000;

ago           = GO47to15(1,1:68)';
agoag         = GO47to15(3,1:68)';
agomi         = GO47to15(6,1:68)';
agout         = GO47to15(10,1:68)';
agoco         = GO47to15(11,1:68)';
agomf         = GO47to15(12,1:68)';
agomfce       = GO47to15(19,1:68)';
agomftemv     = GO47to15(21,1:68)';
agomfteot     = GO47to15(22,1:68)';
agomfch       = GO47to15(32,1:68)';
agowt         = GO47to15(34,1:68)';
agort         = GO47to15(35,1:68)';
agotw         = GO47to15(40,1:68)';
agoin         = GO47to15(49,1:68)';
agoinpu       = GO47to15(50,1:68)';
agoinmp       = GO47to15(51,1:68)';
tem           = GO47to15(52,17)/GO47to15(49,17); 
agoinbt       = [GO47to15(49,1:17)*tem,GO47to15(52,18:68)]';
agoindp       = GO47to15(53,1:68)';
agofr         = GO47to15(54,1:68)';
agofrfi       = GO47to15(55,1:68)';
agopb         = GO47to15(65,1:68)';
agopbpt       = GO47to15(66,1:68)';
agopbptls     = GO47to15(67,1:68)';
agopbptcs     = GO47to15(68,1:68)';
agopbptop     = GO47to15(69,1:68)';
agopbmc       = GO47to15(70,1:68)';
agoeh         = GO47to15(74,1:68)';
agoae         = GO47to15(82,1:68)';
agoos         = GO47to15(89,1:68)';
agogv         = GO47to15(90,1:68)';
agoict        = GO47to15(101,1:68)';

ago15         = [agoag,agomi,agout,agoco,agomf,agowt,agort,agotw, ...
                 agoin,agofr,agopb,agoeh,agoae,agoos,agogv];
ago2          = sum(ago15')';
tgoa          = [1947:2014]';

%
% Impute GO estimates for several intangible-intensive activities
%   (1) Advertising (incl. output in NAICS 54 and 51)
%   (2) R&D (incl. output in NAICS 51, 54 and 31G)
%   (3) Software (incl. output in NAICS 51 and elsewhere)

load ./rawdata/IRS/adexp/irsadvded.dat
tem           = irsadvded(66:-1:1,1:3)/10^6;
tem           = tem(:,2)+[tem(1:42,2)*tem(42,3)/tem(42,2);tem(43:66,3)];
agopbptad     = [tem;tem(66)*agopb(67:68)/agopb(66)];

load ./rawdata/NIPA/FA/table2_7A.dat
load ./rawdata/NIPA/FA/table7_5AA.dat
load ./rawdata/NIPA/FA/table7_5BA.dat

agoinpusf     = [table7_5AA(18,47:96),table7_5BA(20,1:17),NaN]'+ ...
                [table2_7A(79,47:113),table2_7A(79,113)* ...
                 GO403intan(70,18)/GO403intan(70,17)]';
ii            = 1:15;
agoinpusf(ii) = linspace(.5*agoinpusf(15),agoinpusf(15),15)';
agoinpusf(68) = agoinpusf(67)^2/agoinpusf(66);
agopbptrd     = [table7_5AA(19,47:96),table7_5BA(21,1:17),NaN]'+ ...
                [table2_7A(83,47:113),table2_7A(83,113)* ...
                 GO403intan(103,18)/GO403intan(103,17)]';
agopbptrd(68) = agopbptrd(67)^2/agopbptrd(66);

%
% Check imputations with recent data
%

agoinpusfchk  = [NaN*ones(1,50),GO403intan(70,:)]';
agoinbtchk    = [NaN*ones(1,50),sum(GO403intan(73:77,:))]';
agopbptchk    = [NaN*ones(1,50),sum(GO403intan(94:107,:))]';
agopbptlschk  = [NaN*ones(1,50),GO403intan(94,:)]';
agopbptabchk  = [NaN*ones(1,50),GO403intan(98,:)]';
agopbpteachk  = [NaN*ones(1,50),GO403intan(99,:)]';
agopbptcschk  = [NaN*ones(1,50),sum(GO403intan(95:97,:))]';
agopbptmtchk  = [NaN*ones(1,50),sum(GO403intan(101:102,:))]';
agopbptrdchk  = [NaN*ones(1,50),GO403intan(103,:)]';
agopbptadchk  = [NaN*ones(1,50),GO403intan(104,:)]';
agopbptopchk  = [NaN*ones(1,50),sum(GO403intan(105:107,:))]';
agopbmcchk    = [NaN*ones(1,50),GO403intan(108,:)]';

%
% Quaterly estimates start in 2005.1
%

qgo           = [NaN*ones(1,232),GO05to15Q(1,1:40)]';
qgoag         = [NaN*ones(1,232),GO05to15Q(3,1:40)]';
qgomi         = [NaN*ones(1,232),GO05to15Q(4,1:40)]';
qgout         = [NaN*ones(1,232),GO05to15Q(5,1:40)]';
qgoco         = [NaN*ones(1,232),GO05to15Q(6,1:40)]';
qgomf         = [NaN*ones(1,232),GO05to15Q(7,1:40)]';
qgomfce       = [NaN*ones(1,260),GO12to15Qdetail(19,1:12)]';
qgomfee       = [NaN*ones(1,260),GO12to15Qdetail(20,1:12)]';
qgomftemv     = [NaN*ones(1,260),GO12to15Qdetail(21,1:12)]';
qgomfteot     = [NaN*ones(1,260),GO12to15Qdetail(22,1:12)]';
qgomfch       = [NaN*ones(1,260),GO12to15Qdetail(32,1:12)]';
qgowt         = [NaN*ones(1,232),GO05to15Q(10,1:40)]';
qgort         = [NaN*ones(1,232),GO05to15Q(11,1:40)]';
qgotw         = [NaN*ones(1,232),GO05to15Q(12,1:40)]';
qgoin         = [NaN*ones(1,232),GO05to15Q(13,1:40)]';
qgoinpu       = [NaN*ones(1,260),GO12to15Qdetail(50,1:12)]';
qgoinmp       = [NaN*ones(1,260),GO12to15Qdetail(51,1:12)]';
qgoinbt       = [NaN*ones(1,260),GO12to15Qdetail(52,1:12)]';
qgoindp       = [NaN*ones(1,260),GO12to15Qdetail(53,1:12)]';
qgofr         = [NaN*ones(1,232),GO05to15Q(14,1:40)]';
qgofrfi       = [NaN*ones(1,232),GO05to15Q(15,1:40)]';
qgofrre       = [NaN*ones(1,232),GO05to15Q(16,1:40)]';
qgopb         = [NaN*ones(1,232),GO05to15Q(17,1:40)]';
qgopbpt       = [NaN*ones(1,232),GO05to15Q(18,1:40)]';
qgopbptls     = [NaN*ones(1,260),GO12to15Qdetail(67,1:12)]';
qgopbptcs     = [NaN*ones(1,260),GO12to15Qdetail(68,1:12)]';
qgopbptop     = [NaN*ones(1,260),GO12to15Qdetail(69,1:12)]';
qgopbmc       = [NaN*ones(1,232),GO05to15Q(19,1:40)]';
qgoeh         = [NaN*ones(1,232),GO05to15Q(21,1:40)]';
qgoae         = [NaN*ones(1,232),GO05to15Q(24,1:40)]';
qgoos         = [NaN*ones(1,232),GO05to15Q(27,1:40)]';
qgogv         = [NaN*ones(1,232),GO05to15Q(28,1:40)]';
qgo15         = [qgoag,qgomi,qgout,qgoco,qgomf,qgowt,qgort,qgotw,qgoin, ...
                 qgofr,qgopb,qgoeh,qgoae,qgoos,qgogv];
tgoq          = ptime(1947.1,272,4);


logago        = [log(ago15./(defa*ones(1,15))*100),log(agomfch./defa*100), ...
                 log(agoinbt./defa*100),log(agoinpusf./defa*100), ...
                 log(agopbptad./defa*100),log(agomfce./defa*100), ...
                 log(agomftemv./defa*100),log(ago./defa*100)];
tlogago       = [hptrend(logago,100)];

%
% Apply Kalman smoothing routine for choice of business receipt data:
%
%   (0) Total industry
%   (1) Agriculture
%   (2) Mining
%   (3) Utilities
%   (4) Construction
%   (5) Manufacturing
%   (6) Wholesale trade
%   (7) Retail trade
%   (8) Transportation and warehousing
%   (9) Information
%  (10) Finance, insurance, real estate, rental and leasing
%  (11) Professional and business services
%  (12) Educational and health services
%  (13) Arts and entertainment
%  (14) Other services
%  (15) Government
%  (16) Chemical manufacturing
%  (17) Broadcasting and telecommunications
%  (18) Software
%  (19) Advertising 
%  (20) Computer and electronic products
%  (21) Motor vehicles
%
% with all cases covering 1948:1 to 2014:4.
%
%

iest = input('Are quarterly industry data needed? (1=yes, 0=no)' );
if iest==1;
  disp('Annual data available:')
  disp('      1. IRS business receipts ')
  disp('      2. BEA Gross output')
  disp('      3. BEA Value added')
  idata = input(' ? ');
  disp(' ')
  disp('Quarterly data to be computed for which industry:')
  disp('      0. Total industry ')
  disp('      1. Agriculture ')
  disp('      2. Mining ')
  disp('      3. Utilities ')
  disp('      4. Construction ')
  disp('      5. Manufacturing ')
  disp('      6. Wholesale trade ')
  disp('      7. Retail trade ')
  disp('      8. Transportation and warehousing ')
  disp('      9. Information ')
  disp('     10. Finance, information, real estate, rental and leasing ')
  disp('     11. Professional and business services ')
  disp('     12. Educational and health services')
  disp('     13. Arts and entertainment ')
  disp('     14. Other services ')
  disp('     15. Government')
  disp('     16. Chemical manufacturing ')
  disp('     17. Broadcasting and telecommunications')
  disp('     18. Software')
  disp('     19. Advertising ')
  disp('     20. Computer and electronic products ')
  disp('     21. Motor vehicles')
  icase = input('  ?  ');

  if icase==0;
    qni0        = qni(1:268)./defq(5:272);
    qem0        = qem;
    abr0        = [abr;
                   abr(66)*ani(66)/ani(65);
                   abr(66)*ani(67)/ani(65)]./defa;
    ago0        = ago./defa;
    qgo0        = qgo./defq;
  end;
  if icase==1;
    qni0        = qniag(1:268)./defq(5:272);
    qem0        = qemag;
    abr0        = [abrag;
                   abrag(66)*aniag(66)/aniag(65);
                   abrag(66)*aniag(67)/aniag(65)]./defa;
    ago0        = agoag./defa;
    qgo0        = qgoag./defq;
  end;
  if icase==2;
    qni0        = qnimi(1:268)./defq(5:272);
    qem0        = qemmi;
    abr0        = [abrmi;
                   abrmi(66)*animi(66)/animi(65);
                   abrmi(66)*animi(67)/animi(65)]./defa;
    ago0        = agomi./defa;
    qgo0        = qgomi./defq;
  end;
  if icase==3;
    qni0        = qniut(1:268)./defq(5:272);
    qem0        = qemut;
    abr0        = [abrut;
                   abrut(66)*aniut(66)/aniut(65);
                   abrut(66)*aniut(67)/aniut(65)]./defa;
    ago0        = agout./defa;
    qgo0        = qgout./defq;
  end;
  if icase==4;
    qni0        = qnico(1:268)./defq(5:272);
    qem0        = qemco;
    abr0        = [abrco;
                   abrco(66)*anico(66)/anico(65);
                   abrco(66)*anico(67)/anico(65)]./defa;
    ago0        = agoco./defa;
    qgo0        = qgoco./defq;
  end;
  if icase==5;
    qni0        = qnimf(1:268)./defq(5:272);
    qem0        = qemmf;
    abr0        = [abrmf;
                   abrmf(66)*animf(66)/animf(65);
                   abrmf(66)*animf(67)/animf(65)]./defa;
    ago0        = agomf./defa;
    qgo0        = qgomf./defq;
  end;
  if icase==6;
    qni0        = qniwt(1:268)./defq(5:272);
    qem0        = qemwt;
    abr0        = [abrwt;
                   abrwt(66)*aniwt(66)/aniwt(65);
                   abrwt(66)*aniwt(67)/aniwt(65)]./defa;
    ago0        = agowt./defa;
    qgo0        = qgowt./defq;
  end;

  if icase==7;
    qni0        = qnirt(1:268)./defq(5:272);
    qem0        = qemrt;
    abr0        = [abrrt;
                   abrrt(66)*anirt(66)/anirt(65);
                   abrrt(66)*anirt(67)/anirt(65)]./defa;
    ago0        = agort./defa;
    qgo0        = qgort./defq;
  end;
  if icase==8;
    qni0        = qnitw(1:268)./defq(5:272);
    qem0        = qemtw;
    abr0        = [abrtw;
                   abrtw(66)*anitw(66)/anitw(65);
                   abrtw(66)*anitw(67)/anitw(65)]./defa;
    ago0        = agotw./defa;
    qgo0        = qgotw./defq;
  end;
  if icase==9;
    qni0        = qniin(1:268)./defq(5:272);
    qem0        = qemin;
    abr0        = [abrin;
                   abrin(66)*aniin(66)/aniin(65);
                   abrin(66)*aniin(67)/aniin(65)]./defa;
    ago0        = agoin./defa;
    qgo0        = qgoin./defq;
  end;
  if icase==10;
    qni0        = qnifr(1:268)./defq(5:272);
    qem0        = qemfr;
    abr0        = [abrfr;
                   abrfr(66)*anifr(66)/anifr(65);
                   abrfr(66)*anifr(67)/anifr(65)]./defa;
    ago0        = agofr./defa;
    qgo0        = qgofr./defq;
  end;
  if icase==11;
    qni0        = qnipb(1:268)./defq(5:272);
    qem0        = qempb;
    abr0        = [abrpb;
                   abrpb(66)*anipb(66)/anipb(65);
                   abrpb(66)*anipb(67)/anipb(65)]./defa;
    ago0        = agopb./defa;
    qgo0        = qgopb./defq;
  end;
  if icase==12;
    qni0        = qnieh(1:268)./defq(5:272);
    qem0        = qemeh;
    abr0        = [abreh;
                   abreh(66)*anieh(66)/anieh(65);
                   abreh(66)*anieh(67)/anieh(65)]./defa;
    ago0        = agoeh./defa;
    qgo0        = qgoeh./defq;
  end;
  if icase==13;
    qni0        = qniae(1:268)./defq(5:272);
    qem0        = qemae;
    abr0        = [abrae;
                   abrae(66)*aniae(66)/aniae(65);
                   abrae(66)*aniae(67)/aniae(65)]./defa;
    ago0        = agoae./defa;
    qgo0        = qgoae./defq;
  end;
  if icase==14;
    qni0        = qnios(1:268)./defq(5:272);
    qem0        = qemos;
    abr0        = [abros;
                   abros(66)*anios(66)/anios(65);
                   abros(66)*anios(67)/anios(65)]./defa;
    ago0        = agoos./defa;
    qgo0        = qgoos./defq;
  end;
  if icase==15;
    qni0        = qnigv(1:268)./defq(5:272);
    qem0        = qemgv;
    ago0        = agogv./defa;
    qgo0        = qgogv./defq;
    abr0        = NaN*ago0;
  end;


  if icase==16;
    qni0        = qnimf(1:268)./defq(5:272);
    qem0        = qemmfch;
    abr0        = [abrmfch;
                   abrmfch(66)*animf(66)/animf(65);
                   abrmfch(66)*animf(67)/animf(65)]./defa;
    ago0        = agomfch./defa;
  end;

  if icase==17;
    qni0        = qniin(1:268)./defq(5:272);
    
    qem0        = qeminbt;
    abr0        = [abrinbt;
                   abrinbt(66)*aniin(66)/aniin(65);
                   abrinbt(66)*aniin(67)/aniin(65)]./defa;
    ago0        = agoinbt./defa;
  end;

  if icase==18;
    qni0        = qniin(1:268)./defq(5:272);
    qem0        = qeminpusf;
    ago0        = agoinpusf./defa;
    abr0        = NaN*ago0;
  end;

  if icase==19;
    qni0        = qnipb(1:268)./defq(5:272);
    % 
    % Note: found strange oscillations using qnipb so used
    %       two relevant employment series
    % 
    qni0        = qeminbtrt;
    ii          = 1:44;
    qni0(ii)    = qemin(1:44)*qni0(45)/qemin(45);
    qni0        = qni0(5:272);

    qem0        = qempbptad;
    ii          = 1:44;
    qem0(ii)    = qempb(1:44)*qem0(45)/qempb(45);

    abr0        = [abrpbptad;
                   abrpbptad(66)*anipb(66)/anipb(65);
                   abrpbptad(66)*anipb(67)/anipb(65)]./defa;
    ago0        = agopbptad./defa;
  end;
  if icase==20;
    qni0        = qnimf(1:268)./defq(5:272);
    qem0        = qemmfce;
    ago0        = agomfce./defa;
    abr0        = NaN*ago0;
  end;
  if icase==21;
    qni0        = qnimf(1:268)./defq(5:272);
    qem0        = qemmftemv;
    ago0        = agomftemv./defa;
    abr0        = NaN*ago0;
  end;

  tqni0         = hptrend(log(qni0),1600);
  tqem0         = hptrend(log(qem0),1600);
  qem0          = qem0(5:272);
  tqem0         = tqem0(5:272);
  tabr0         = hptrend(log(abr0),100);
  abr0          = abr0(2:68);
  tabr0         = tabr0(2:68);
  tago0         = hptrend(log(ago0),100);
  ago0          = ago0(2:68);
  tago0         = tago0(2:68);

  ii            = 4:4:268;
  n             = 3;
  Y             = zeros(268,n);
  Y(:,1:2)      = [log(qni0)-tqni0,log(qem0)-tqem0];

  %if icase==18;
  %  n           = 2; 
  %  Y           = zeros(268,n);
  %  Y(:,1)      = log(qem0)-tqem0;
  %end; 

  if idata==1;
    Y(ii,n)     = log(abr0)-tabr0;
    qYbr        = ann2qtr([[1948:2014]',Y(4:4:268,n)],4,1);
    qYbr        = qYbr(:,2);
    qY          = [Y(:,1:n-1),qYbr];
  else;
    Y(ii,n)     = log(ago0)-tago0;
    qYgo        = ann2qtr([[1948:2014]',Y(4:4:268,n)],4,1);
    qYgo        = qYgo(:,2);
    qY          = [Y(:,1:n-1),qYgo];
  end;

  sumyx         = zeros(n,4*n+1);
  sumxx         = zeros(4*n+1,4*n+1);
  for t=5:268;
    Ylag        = [qY(t-1,:)';qY(t-2,:)';qY(t-3,:)';qY(t-4,:)';1];
    sumyx       = sumyx+qY(t,:)'*Ylag';
    sumxx       = sumxx+Ylag*Ylag';
  end;
  A             = sumyx*inv(sumxx); 
  Sig           = zeros(n);
  for t=5:268;
    Ylag        = [qY(t-1,:)';qY(t-2,:)';qY(t-3,:)';qY(t-4,:)';1];
    eps         = qY(t,:)'-A*Ylag;
    Sig         = Sig+eps*eps'/264;
  end;
  A             = A(:,1:4*n);
  B             = chol(Sig)';
  Theta0        = [A(:);B(:)];


  if (iest==1);
 
    [Theta1,f,g]  = uncmin(Theta0,'lestqtr',Y)
    if abs(sum(imag(Theta1)))>0;
      disp('Warning: imaginary numbers found in Theta')
      Y2          = estqtr(real(Theta1),Y);
    else
      Y2          = estqtr(Theta1,Y);
    end;
    if (idata==1)
      qtrend    = ann2qtr([[1948:2014]',tabr0],4,1);
    else
      qtrend    = ann2qtr([[1948:2014]',tago0],4,1);
    end;
    qtrend      = qtrend(:,2);
    qgopred     = exp(Y2(:,n)+qtrend);
    if (idata==2 & icase<16);
      %
      % Use quarterly gross output figures from BEA
      %
      qgoact    = qgo0(5:272);
      qgopred   = [qgopred(1:229)*qgoact(229)/qgopred(229);qgoact(230:268)];
    end;
    Y2n         = Y2;
    Y2n(:,n)    = log(qgopred)-hptrend(log(qgopred),1600);
  end;

end;

