%INPUTS   this file loads in files with data on aggregate FDI flows, GDPs, 
%         and populations and files with bilateral inward and outward FDI
%         flows. To include more countries in the analysis, users need to 
%         update the files being loaded.  (See ./data/README for more details.)
%

%         Ellen McGrattan, 10-23-12
%         Revised, 4-23-13, ERM


clear
%============================================================================
% LOAD AGGREGATE DATA 
%
%   * GDP and per capita GDP, current USD (NGDPD/NGDPDPC)
%   * GDP and per capita GDP, constant USD (RGDPD/RGDPDPC)
%   * GDP and per capita GDP, current PPP (NGDPP/NGDPPPC)
%   * GDP and per capita GDP, constant PPP (RGDPP/RGDPPC)
%   * Total population (POP)
%   * Population ages 15-64 (POP1564)
%   * Inward FDI flows (FDII)
%   * Corporate tax rates

%
%   NOTE: The aggregate data files loaded into this code have a format 
%         with years across columns and countries across rows:
%
%                 1. Years (1990-2010 available)
%
%                 2. US                                
%                                                      ____
%                 3. Austria                               |
%                 4. Belgium+Luxembourg                    |
%                 5. Denmark                               |  EU-16
%                 6. Finland                               |
%                 7. France                                |
%                 8. Germany                               |
%                 9. Greece                                |
%                10. Ireland                               |
%                11. Italy                                 |
%                12. Netherlands                           |
%                13. Norway                                |
%                14. Portugal                              |
%                15. Spain                                 |
%                16. Sweden                                |
%                17. Switzerland                           |
%                18. United Kingdom                        |
%                                                      ----
%                19. Japan
%                                                      ____
%                20. Korea                                 | OECD that
%                21. Australia                             | are in ROW
%                22. Canada                            ----
%                                                      ____
%                23. China                                 | 
%                24. Hong Kong                             | CHT of CHTM
%                25. Taiwan                            ----
%                                                      ____
%                26. India                                 |
%                27. Russia                                | BRI
%                28. Brazil                            ----
%                                                      ____
%                29. Bahamas                               | BBB
%                30. Bermuda                               |
%                31. Barbados                          ----
%                                                      ____
%                32. Cayman Islands                        | 
%                33. BVI                                   | UKI
%                34. UK Islands (includes Cayman and BVI   |
%                                if not in rows 32-33) ----
%
%                35. Bahrain                               | Part of ROW
%
%                36. Macao                                 | M of CHTM
%                                                      ____
%                37. Indonesia                             |
%                38. Malaysia                              |
%                39. Mauritius                             |
%                40. Philippines                           | Rest of
%                41. Samoa                                 | ROW
%                42. Saudi Arabia                          |
%                43. Singapore                             |
%                44. Thailand                          ----
%
%              To add more countries, raw data files need to be edited.
%              (See main README in ./data for details.)
%
%============================================================================

t0             = 1990;
t1             = 2010;
t              = t0:t1;
n              = 44;
label          = ['Year';' USA';' AUT';' BLX';' DNK';' FIN';' FRA';' DEU';
                  ' GRC';' IRL';' ITA';' NLD';' NOR';' PRT';' ESP';' SWE';
                  ' CHE';' GBR';' JPN';' KOR';' AUS';' CAN';' CHN';' HKG';
                  ' TWN';' IND';' RUS';' BRA';' BHS';' BMU';' BRB';' CYM';
                  ' BVI';' UKI';' BHR';' MAC';' IDN';' MYS';' MUS';' PHL';
                  ' WSM';' SAU';' SGP';' THA'];
load rawdata/unsd/unsd.out
un_ngdpd       = unsd(1:n,:);
un_rgdpd       = unsd(n+1:2*n,:);
un_ngdpdpc     = unsd(2*n+1:3*n,:);
un_xrate       = unsd(3*n+1:4*n,:);
un_poptot      = unsd(4*n+1:5*n,:);

load rawdata/wdi/wdi.out
wdi_ngdpd      = wdi(1:n,:);
wdi_ngdpp      = wdi(n+1:2*n,:);
wdi_rgdpd      = wdi(2*n+1:3*n,:);
wdi_rgdpp      = wdi(3*n+1:4*n,:);
wdi_ngdpdpc    = wdi(4*n+1:5*n,:);
wdi_ngdpppc    = wdi(5*n+1:6*n,:);
wdi_rgdpdpc    = wdi(6*n+1:7*n,:);
wdi_rgdpppc    = wdi(7*n+1:8*n,:);
wdi_poptot     = wdi(8*n+1:9*n,:);
wdi_pop1564    = wdi(9*n+1:10*n,:);
wdi_cabal      = wdi(10*n+1:11*n,:);
wdi_expgs      = wdi(11*n+1:12*n,:);
wdi_impgs      = wdi(12*n+1:13*n,:);
wdi_bopincr    = wdi(13*n+1:14*n,:);
wdi_bopincp    = wdi(14*n+1:15*n,:);
wdi_trannet    = wdi(15*n+1:16*n,:);
wdi_capnet     = wdi(16*n+1:17*n,:);
wdi_errnet     = wdi(17*n+1:18*n,:);
wdi_fdiin      = wdi(18*n+1:19*n,:);
wdi_fdinet     = wdi(19*n+1:20*n,:);
wdi_pinet      = wdi(20*n+1:21*n,:);
wdi_ictexp     = wdi(21*n+1:22*n,:);
wdi_royrec     = wdi(22*n+1:23*n,:);
wdi_roypay     = wdi(23*n+1:24*n,:);

load rawdata/weo/weo.out
weo_ngdpd      = [weo(1,:);weo(2:n,:)*10^9];
weo_ngdpp      = [weo(n+1,:);weo(n+2:2*n,:)*10^9];
weo_ngdpdpc    = [weo(2*n+1,:);weo(2*n+2:3*n,:)];
weo_ngdpppc    = [weo(3*n+1,:);weo(3*n+2:4*n,:)];
weo_poptot     = [weo(4*n+1,:);weo(4*n+2:5*n,:)*10^6];

load rawdata/ted/ted.out
ted_rgdpe     = ted(1:n,:);
ted_rgdpg     = ted(n+1:2*n,:);
ted_rgdpepc   = ted(2*n+1:3*n,:);
ted_rgdpgpc   = ted(3*n+1:4*n,:);
ted_poptot    = ted(4*n+1:5*n,:);

load rawdata/imf/imf.out
imf_fdii      = [imf(1,:);imf(2:n,:)*10^6];

load rawdata/unctad/unctad.out
load rawdata/unctad/unctado.out
un_fdii       = [unctad(1,:);unctad(2:n,:)*10^6];
un_fdio       = [unctado(1,:);unctado(2:n,:)*10^6];

load rawdata/oecd/taxes/corptax/oecd.out
load rawdata/kpmg/taxes/kpmg.out
taup          = [oecd(1:22,1:3),kpmg(1:22,4:length(t));
                 kpmg(23:44,8)*ones(1,7),kpmg(23:44,8:length(t))];

%
% WDI is default for GDPs and population unless data are missing
%
ngdpd         = wdi_ngdpd;
ngdpp         = wdi_ngdpp;
rgdpd         = wdi_rgdpd;
rgdpp         = wdi_rgdpp;
ngdpdpc       = wdi_ngdpdpc;
ngdpppc       = wdi_ngdpppc;
rgdpdpc       = wdi_rgdpdpc;
rgdpppc       = wdi_rgdpppc;
poptot        = wdi_poptot;

%
% Fill in missing data for Ireland
%
t00           = find(t==2000);
t05           = find(t==2005);
rgdpd(10,:)   = un_rgdpd(10,:) * wdi_rgdpd(10,t00)/un_rgdpd(10,t00);
rgdpp(10,:)   = ted_rgdpg(10,:) * wdi_rgdpp(10,t05)/ted_rgdpg(10,t05);
rgdpdpc(10,:) = rgdpd(10,:)./poptot(10,:);
rgdpppc(10,:) = rgdpp(10,:)./poptot(10,:);

%
% Fill in missing data for Taiwan
%
ngdpd(25,:)   = weo_ngdpd(25,:);
ngdpp(25,:)   = weo_ngdpp(25,:);
ngdpdpc(25,:) = weo_ngdpdpc(25,:);
ngdpppc(25,:) = weo_ngdpppc(25,:);
poptot(25,:)  = weo_poptot(25,:);
rgdpd(25,:)   = ted_rgdpg(25,:) *ngdpd(25,t00)/ted_rgdpg(25,t00);
rgdpp(25,:)   = ted_rgdpg(25,:) *ngdpp(25,t05)/ted_rgdpg(25,t05);
rgdpdpc(25,:) = rgdpd(25,:)./poptot(25,:);
rgdpppc(25,:) = rgdpp(25,:)./poptot(25,:);

%
% Fill in missing data for Bermuda 
%
ngdpp(30,:)   = ngdpd(30,:);
rgdpp(30,:)   = rgdpd(30,:) *ngdpd(30,t05)/rgdpd(30,t05);
ngdpppc(30,:) = ngdpp(30,:)./poptot(30,:);
rgdpppc(30,:) = rgdpp(30,:)./poptot(30,:);

%
% Fill in missing data for Barbados
%
t10           = find(t==2010);
ngdpp(31,t10) = weo_ngdpp(31,t10)*ngdpp(31,t10-1)/weo_ngdpp(31,t10-1);
rgdpd(31,t10) = un_rgdpd(31,t10)*rgdpd(31,t10-1)/un_rgdpd(31,t10-1);
rgdpp(31,t10) = rgdpp(31,t10-1)*rgdpd(31,t10)/rgdpd(31,t10-1);
ngdpppc(31,:) = ngdpp(31,:)./poptot(31,:);
rgdpdpc(31,:) = rgdpd(31,:)./poptot(31,:);
rgdpppc(31,:) = rgdpp(31,:)./poptot(31,:);


%
% Fill in missing data for UK islands and Macao
%
ii            = [32:34,36]; 
nt            = length(t);
poptot(ii,:)  = un_poptot(ii,:);
ngdpd(ii,:)   = un_ngdpd(ii,:);
ngdpp(ii,:)   = un_ngdpd(ii,:);   
ngdpdpc(ii,:) = un_ngdpdpc(ii,:);
ngdpppc(ii,:) = un_ngdpdpc(ii,:);
rgdpd(ii,:)   = un_rgdpd(ii,:).*((ngdpd(ii,t00)./un_rgdpd(ii,t00))*ones(1,nt));
rgdpp(ii,:)   = rgdpd(ii,:).*((ngdpp(ii,t05)./rgdpd(ii,t05))*ones(1,nt));
rgdpdpc(ii,:) = rgdpd(ii,:)./poptot(ii,:);
rgdpppc(ii,:) = rgdpp(ii,:)./poptot(ii,:);

%
% UNCTAD: default for FDI flows
%
fdii          = un_fdii;

%
% USD deflator
%
pgdpd         = ngdpd./rgdpd;


%============================================================================
% USER-DEFINED COUNTRY GROUPINGS AND SAMPLES
% 
%   NOTE: The default choices, which can be changed, are given by:
%
%            (1) Country groupings, where the default is:
%                Region 1 = US
%                       2 = EU-16 (see list below) 
%                       3 = Japan+Korea (JK)
%                       4 = China+Provinces (CHTM)
%                       5 = Brazil+Russia+India (BRI)
%                       6 = ROW significant for China (see list below)
%                       7 = Bahamas, Bermuda, Barbados (BBB)
%                       8 = UK islands (UKI)
%                       9 = Regions 1:6
%                      10 = Regions 1:8
%
%            (2) Samples, where the default is: 
%                Period 1 = 1990-1990
%                       2 = 1991-1995
%                       3 = 1996-2000
%                       4 = 2001-2005
%                       5 = 2006-2010
%                       6 = 2001-2010 (proxy for long-run)
%

nr            = 10;
region1       = 2;
region2       = 3:18;
region3       = 19:20;
region4       = [23:25,36];
region5       = 26:28;
region6       = [21:22,35,37:44];
region7       = 29:31;
region8       = 32:34;
region9       = [region1,region2,region3,region4,region5,region6];
region10      = [region9,region7,region8];

np            = 6;
period1       = [1990,1990];
period2       = [1991,1995];
period3       = [1996,2000];
period4       = [2001,2005];
period5       = [2006,2010];
period6       = [2001,2010];


%=============================================================================
% STATISTICS OF INTEREST 
%
%   * Bilateral FDI inward flows for US, EU-16, Japan+Korea, China & provinces
%   * Inward FDI relative to GDP (USD)
%   * Total population relative to US
%   * Real per capita GDP (PPP) relative to US
%   * Real per capita FDI (USD) relative to a 2% trend
%   * Real per capita GDP (USD) relative to a 2% trend
%   * Corporate tax rates, population weighted
%

load rawdata/bilatflows/other/macau/bilat_fdii.out
tt            = bilat_fdii(1,:);
i             = find(tt==t0);
j             = find(tt==t1);
ximac         = [bilat_fdii(1,i:j);
                 bilat_fdii(2:45,i:j)*10^6;
                 un_fdii(36,:)];
ximac(45,1:11)= ximac(46,1:11);

load rawdata/bilatflows/other/macau/bilat_fdio.out
tt            = bilat_fdio(1,:);
i             = find(tt==t0);
j             = find(tt==t1);
xomac         = [bilat_fdio(1,i:j);
                 bilat_fdio(2:45,i:j)*10^6;
                 un_fdio(36,:)];
xomac(45,1:11)= xomac(46,1:11);


sum1          = 0.;
sum2          = 0.;
sum3          = 0.;
sum4          = 0.;
sum5          = 0.;
sum6          = 0.;
sum4o         = 0.;
for i=t0:t1;
  s           = int2str(i);
  eval(['load rawdata/bilatflows/inflows/year',s,'.txt'])
  eval(['year',s,'i = year',s,';'])
  eval(['xi = year',int2str(i),';'])
  eval(['load rawdata/bilatflows/outflows/year',s,'.txt'])
  eval(['year',s,'o = year',s,';'])
  eval(['xo = year',int2str(i),';'])

  xi          = [NaN*ones(46,1),[i*ones(1,23);xi*10^6]];
  xo          = [NaN*ones(46,1),[i*ones(1,23);xo*10^6]];
  j           = i-t0+1;

  %
  % Use outward flows for Taiwan
  %
  xitwn(:,j)  = [i;xo(25,2:24)';NaN*ones(20,1);
                 un_fdii(25,j);un_fdii(25,j)];
  xotwn(:,j)  = [i;xi(25,2:24)';NaN*ones(20,1);
                 un_fdio(25,j);un_fdio(25,j)];

  %
  % Data for Hong Kong are unavailable 1990-1997. We have 
  % total flows from Unctad and bilateral stocks reported by
  % Hong Kong in 1998. We use these stocks to estimate bilateral
  % flows prior to 1998.
  %
  if (i<1998);
    xi(2,24)  = 108.3/1627.2*un_fdii(24,j);
    xi(12,24) = 119.9/1627.2*un_fdii(24,j);
    xi(18,24) = 149.3/1627.2*un_fdii(24,j);
    xi(19,24) = 105.5/1627.2*un_fdii(24,j);
    xi(23,24) = 205.9/1627.2*un_fdii(24,j);
    xi(30,24) = 198.1/1627.2*un_fdii(24,j);
    xi(32,24) = 77.8/1627.2*un_fdii(24,j);
    xi(33,24) = 514.9/1627.2*un_fdii(24,j);
    xi(43,24) = 37.5/1627.2*un_fdii(24,j);
    xi(45,24) = un_fdii(24,j);
    xi(46,24) = un_fdii(24,j);

    xo(2,24)  = 12.9/1668.8*un_fdio(24,j);
    xo(18,24) = 92.0/1668.8*un_fdio(24,j);
    xo(21,24) = 19.8/1668.8*un_fdio(24,j);
    xo(23,24) = 534.3/1668.8*un_fdio(24,j);
    xo(30,24) = 88.1/1668.8*un_fdio(24,j);
    xo(32,24) = 89.5/1668.8*un_fdio(24,j);
    xo(33,24) = 692.6/1668.8*un_fdio(24,j);
    xo(44,24) = 7.6/1668.8*un_fdio(24,j);
    xo(45,24) = un_fdio(24,j);
    xo(46,24) = un_fdio(24,j);

  end;

  %
  % Outward FDI flows for China are unavailable 1990-2002. We have 
  % total flows from Unctad and bilateral stocks reported by
  % for China in 2003. We use these stocks to estimate bilateral
  % flows prior to 2003.
  %
  if (i<2003);

    xo(2,23)  = 502.32/33222.22*un_fdio(23,j);
    xo(3,23)  = 0.70/33222.22*un_fdio(23,j);
    xo(4,23)  = 0.41/33222.22*un_fdio(23,j);
    xo(5,23)  = 74.43/33222.22*un_fdio(23,j);
    xo(7,23)  = 13.12/33222.22*un_fdio(23,j);
    xo(8,23)  = 83.61/33222.22*un_fdio(23,j);
    xo(10,23) = 0.24/33222.22*un_fdio(23,j);
    xo(11,23) = 19.18/33222.22*un_fdio(23,j);
    xo(12,23) = 5.90/33222.22*un_fdio(23,j);
    xo(14,23) = 0.20/33222.22*un_fdio(23,j);
    xo(15,23) = 101.81/33222.22*un_fdio(23,j);
    xo(16,23) = 6.07/33222.22*un_fdio(23,j);
    xo(17,23) = 2.21/33222.22*un_fdio(23,j);
    xo(18,23) = 75.15/33222.22*un_fdio(23,j);
    xo(19,23) = 89.31/33222.22*un_fdio(23,j);
    xo(20,23) = 1.17/33222.22*un_fdio(23,j);
    xo(21,23) = 416.49/33222.22*un_fdio(23,j);
    xo(22,23) = 46.8/33222.22*un_fdio(23,j);
    xo(24,23) = 24632.26/33222.22*un_fdio(23,j);
    xo(26,23) = 0.96/33222.22*un_fdio(23,j);
    xo(27,23) = 61.64/33222.22*un_fdio(23,j);
    xo(28,23) = 52.19/33222.22*un_fdio(23,j);
    xo(29,23) = 44.45/33222.22*un_fdio(23,j);
    xo(32,23) = 3690.68/33222.22*un_fdio(23,j);
    xo(33,23) = 532.64/33222.22*un_fdio(23,j);
    xo(35,23) = 0.15/33222.22*un_fdio(23,j);
    xo(36,23) = 444.86/33222.22*un_fdio(23,j);
    xo(37,23) = 54.26/33222.22*un_fdio(23,j);
    xo(38,23) = 100.66/33222.22*un_fdio(23,j);
    xo(39,23) = 12.59/33222.22*un_fdio(23,j);
    xo(40,23) = 8.75/33222.22*un_fdio(23,j);
    xo(41,23) = 0.92/33222.22*un_fdio(23,j);
    xo(42,23) = 0.24/33222.22*un_fdio(23,j);
    xo(43,23) = 164.83/33222.22*un_fdio(23,j);
    xo(44,23) = 150.77/33222.22*un_fdio(23,j);
    xo(45,23) = un_fdio(23,j);
    xo(46,23) = un_fdio(23,j);

  end;

  b1(:,j)     = xi(:,region1);
  b2(:,j)     = sumnan(xi(:,region2)')';
  b3(:,j)     = sumnan(xi(:,region3)')';
  b4(:,j)     = sumnan([xi(:,23:24),xitwn(:,j),ximac(:,j)]')';

  for l=1:4
    lstr       = int2str(l);
    eval(['own',lstr,'=sumnan(b',lstr,'(region',lstr,',j));']);
    for k=1:nr;
      kstr    = int2str(k);
      if (l~=k)
        eval(['b',kstr,'to',lstr,'=sumnan(b',lstr,'(region',kstr,',j));']);

      end;
    end;
  end;

  bwto1       = b1(45,j);
  bwto2       = b2(45,j)-own2;
  bwto3       = b3(45,j)-own3;
  bwto4       = b4(45,j)-own4;
  bwto5       = sumnan(fdii(region5,j));
  bwto6       = sumnan(fdii(region6,j));

  bito1       = imf_fdii(region1,j);
  bito2       = sumnan(imf_fdii(region2,j))-own2;
  bito3       = sumnan(imf_fdii(region3,j))-own3;
  bito4       = sumnan(imf_fdii(region4,j))-own4;
  bito5       = sumnan(imf_fdii(region5,j));
  bito6       = sumnan(imf_fdii(region6,j));

  buto1       = un_fdii(region1,j);
  buto2       = sumnan(un_fdii(region2,j))-own2;
  buto3       = sumnan(un_fdii(region3,j))-own3;
  buto4       = sumnan(un_fdii(region4,j))-own4;
  buto5       = sumnan(un_fdii(region5,j));
  buto6       = sumnan(un_fdii(region6,j));

  bbto1       = wdi_fdiin(region1,j);
  bbto2       = sumnan(wdi_fdiin(region2,j))-own2;
  bbto3       = sumnan(wdi_fdiin(region3,j))-own3;
  bbto4       = sumnan(wdi_fdiin(region4,j))-own4;
  bbto5       = sumnan(wdi_fdiin(region5,j));
  bbto6       = sumnan(wdi_fdiin(region6,j));

  b9to2       = b9to2-own2;
  b9to4       = b9to4-own4;
  b10to2      = b10to2-own2;
  b10to4      = b10to4-own4;

  %
  % Outward from China
  %
  b4o(:,j)    = sumnan([xo(:,23:24),xotwn(:,j),xomac(:,j)]')';
  lstr        = int2str(4);
  own4o       = sumnan(b4o(region4,j));
  for k=1:nr;
    kstr      = int2str(k);
    if (k~=4)
      eval(['b',kstr,'from',lstr,'=sumnan(b4o(region',kstr,',j));']);
    end;
  end;
  bwfrom4     = b4o(45,j)-own4o;

  for k=1:nr
    kstr      = int2str(k);
    eval(['l  = region',int2str(k),';']);
    if (length(l)>1);
      eval(['pop',  kstr,'=sumnan(poptot(l,j))/poptot(2,j)*100;']);
      eval(['taup', kstr,'=sumnan(poptot(l,j).*taup(l,j))/', ...
            'sumnan(poptot(l,j).*(~isnan(taup(l,j))));']);
      eval(['gdppc',kstr,'=(sumnan(rgdpp(l,j))/sumnan(poptot(l,j)))./', ...
                          '(sumnan(rgdpp(2,j))/sumnan(poptot(2,j)))*100;']);

      eval(['pgdpd',kstr,'=sumnan(ngdpd(l,j))/sumnan(rgdpd(l,j)) *', ...
                          'sumnan(rgdpd(l,1))/sumnan(ngdpd(l,1));']);
      eval(['ngdp',kstr,'=sumnan(ngdpd(l,j));']);
      eval(['trend',kstr,'=(pgdpd',kstr,'*sumnan(poptot(l,j))*1.02^(j-1));']);
      eval(['dgdp',kstr,'=sumnan(ngdpd(l,j))/trend',kstr,';']);
    
    else;
      eval(['pop',  kstr,'=poptot(l,j)/poptot(2,j)*100;']);
      eval(['taup', kstr,'=taup(l,j);']);
      eval(['gdppc',kstr,'=(rgdpp(l,j)/poptot(l,j))./', ...
                          '(rgdpp(2,j)/poptot(2,j))*100;'])
      eval(['pgdpd',kstr,'=ngdpd(l,j)/rgdpd(l,j) * rgdpd(l,1)/ngdpd(l,1);']);
      eval(['ngdp',kstr,'=ngdpd(l,j);']);
      eval(['trend',kstr,'=(pgdpd',kstr,'*poptot(l,j)*1.02^(j-1))      ;']);
      eval(['dgdp',kstr,'=ngdpd(l,j)/trend',kstr,';']);
    end;
  end;

  %
  % Inward FDI statistics
  %

  fdigdp1a    = bwto1/ngdp1;
  fdigdp2a    = bwto2/ngdp2;
  fdigdp3a    = bwto3/ngdp3;
  fdigdp4a    = bwto4/ngdp4;
  fdigdp1b    = sumnan([bwto1;-b7to1;-b8to1])/ngdp1;
  fdigdp2b    = sumnan([bwto2;-b7to2;-b8to2])/ngdp2;
  fdigdp3b    = sumnan([bwto3;-b7to3;-b8to3])/ngdp3;
  fdigdp4b    = sumnan([bwto4;-b7to4;-b8to4])/ngdp4;
  fdigdp1c    = sumnan([bwto1;-b7to1])/ngdp1;
  fdigdp2c    = sumnan([bwto2;-b7to2;-b8to2])/ngdp2;
  fdigdp3c    = bwto3/ngdp3;
  fdigdp4c    = sumnan([bwto4;-b8to4])/ngdp4;
  fdigdp5     = sumnan( (un_fdii(region5,j).*poptot(region5,j))./ ...
                ngdpd(region5,j))/sumnan(poptot(region5,j));
  fdigdp6     = sumnan( (un_fdii(region6,j).*poptot(region6,j))./ ...
                ngdpd(region6,j))/sumnan(poptot(region6,j));
  %
  % Outward FDI statistics from CHTM
  %
  ofdigdp4  = sumnan([bwfrom4;-b8from4])/ngdp4;

  %
  % Net portfolio flows (new liabilities at home + neg. new assets abroad)
  %
  pinet1      = sumnan(wdi_pinet(region1,j));
  pinet2      = sumnan(wdi_pinet(region2,j));
  pinet3      = sumnan(wdi_pinet(region3,j));
  pinet4      = sumnan(wdi_pinet(region4,j));
  pinet5      = sumnan(wdi_pinet(region5,j));
  pinet6      = sumnan(wdi_pinet(region6,j));
  pingdp1     = pinet1/ngdp1;
  pingdp2     = pinet2/ngdp2;
  pingdp3     = pinet3/ngdp3;
  pingdp4     = pinet4/ngdp4;
  pingdp5     = pinet5/ngdp5;
  pingdp6     = pinet6/ngdp6;

  %
  % Which FDI/GDP series?
  %
  icase       = 3;
  if (icase==1);
    fdigdp1   = fdigdp1a;
    fdigdp2   = fdigdp2a;
    fdigdp3   = fdigdp3a;
    fdigdp4   = fdigdp4a;
  elseif (icase==2);
    fdigdp1   = fdigdp1b;
    fdigdp2   = fdigdp2b;
    fdigdp3   = fdigdp3b;
    fdigdp4   = fdigdp4b;
  else
    fdigdp1   = fdigdp1c;
    fdigdp2   = fdigdp2c;
    fdigdp3   = fdigdp3c;
    fdigdp4   = fdigdp4c;
  end;
  dfdi1       = fdigdp1*ngdp1/trend1;
  dfdi2       = fdigdp2*ngdp2/trend2;
  dfdi3       = fdigdp3*ngdp3/trend3;
  dfdi4       = fdigdp4*ngdp4/trend4;
  dfdi5       = fdigdp5*ngdp5/trend5;
  dfdi6       = fdigdp6*ngdp6/trend6;
  sum1        = sumnan([sum1;fdigdp1*ngdp1]);
  sum2        = sumnan([sum2;fdigdp2*ngdp2]);
  sum3        = sumnan([sum3;fdigdp3*ngdp3]);
  sum4        = sumnan([sum4;fdigdp4*ngdp4]);
  sum5        = sumnan([sum5;fdigdp5*ngdp5]);
  sum6        = sumnan([sum6;fdigdp6*ngdp6]);
  dfdik1      = sum1 /trend1;
  dfdik2      = sum2 /trend2;
  dfdik3      = sum3 /trend3;
  dfdik4      = sum4 /trend4;
  dfdik5      = sum5 /trend5;
  dfdik6      = sum6 /trend6;

  sum4o       = sumnan([sum4o;ofdigdp4*ngdp4]);
  dfdik4o     = sum4o /trend4;

  disp(' ')
  disp(' ')
  disp(' ')
  disp('==========')
  disp(sprintf('YEAR: %4g',t(j)))
  disp('==========')
  disp('                     Inward FDI To:         | Relative to US:|       Detrended:       | %FDI/GDP      ')
  disp('                US      EU    Jap-Kor   CHTM|  Pop    GDP pc |    FDI  FDI-K   GDP    | w, wo/ BBB+UKI')
  disp('--------------------------------------------|----------------|------------------------|---------------')
  disp('From:                                       |                |                        |')
  disp(sprintf('  US               %9.1f %7.1f %7.1f|%7.1f %7.1f |%7.1f %7.1f %7.1f | %6.3f %6.3f', ...
  [[b1to2,b1to3,b1to4]/10^9,pop1,gdppc1,dfdi1,dfdik1,dgdp1,fdigdp1a*100,fdigdp1b*100]));
  disp(sprintf('  EU-16      %7.1f       %9.1f %7.1f|%7.1f %7.1f |%7.1f %7.1f %7.1f | %6.3f %6.3f', ... 
  [[b2to1,b2to3,b2to4]/10^9,pop2,gdppc2,dfdi2,dfdik2,dgdp2,fdigdp2a*100,fdigdp2b*100]));
  disp(sprintf('  Jap-Kor    %7.1f %7.1f       %9.1f|%7.1f %7.1f |%7.1f %7.1f %7.1f | %6.3f %6.3f', ... 
  [[b3to1,b3to2,b3to4]/10^9,pop3,gdppc3,dfdi3,dfdik3,dgdp3,fdigdp3a*100,fdigdp3b*100]));
  disp(sprintf('  CHTM       %7.1f %7.1f %7.1f        |%7.1f %7.1f |%7.1f %7.1f %7.1f | %6.3f %6.3f', ... 
  [[b4to1,b4to2,b4to3]/10^9,pop4,gdppc4,dfdi4,dfdik4,dgdp4,fdigdp4a*100,fdigdp4b*100]));
  disp(sprintf('  BRI        %7.1f %7.1f %7.1f %7.1f|%7.1f %7.1f |%7.1f %7.1f %7.1f | %6.3f', ... 
  [[b5to1,b5to2,b5to3,b5to4]/10^9,pop5,gdppc5,dfdi5,dfdik5,dgdp5,fdigdp5*100]));
  disp(sprintf('  ROW-12     %7.1f %7.1f %7.1f %7.1f|%7.1f %7.1f |%7.1f %7.1f %7.1f | %6.3f', ... 
  [[b6to1,b6to2,b6to3,b6to4]/10^9,pop6,gdppc6,dfdi6,dfdik6,dgdp6,fdigdp6*100]));
  disp('                                            |                |')
  disp(sprintf('  BBB        %7.1f %7.1f %7.1f %7.1f|%7.1f %7.1f |', ...
  [[b7to1,b7to2,b7to3,b7to4]/10^9,pop7,gdppc7]));
  disp(sprintf('  UKI        %7.1f %7.1f %7.1f %7.1f|%7.1f %7.1f |', ...
  [[b8to1,b8to2,b8to3,b8to4]/10^9,pop8,gdppc8]));
  disp('                                            |      ')
  disp(sprintf('  Residual   %7.1f %7.1f %7.1f %7.1f|      ', ...
  sumnan([bwto1,bwto2,bwto3,bwto4;
        -[b10to1,b10to2,b10to3,b10to4]]/10^9)));
  disp('---------------------------------------------')
  disp(sprintf('  Rows 1-6   %7.1f %7.1f %7.1f %7.1f  ', ...
  [b9to1,b9to2,b9to3,b9to4]/10^9));
  disp(sprintf('  Rows 1-8   %7.1f %7.1f %7.1f %7.1f  ', ...
  [b10to1,b10to2,b10to3,b10to4]/10^9));
  disp('--------------------------------------------')
  disp(sprintf('World ($B)   %7.1f %7.1f %7.1f %7.1f  ', ...
  [bwto1,bwto2,bwto3,bwto4]/10^9))
  disp(sprintf('  Less BBB   %7.1f %7.1f %7.1f %7.1f  ', ...
  sumnan([bwto1,bwto2,bwto3,bwto4;
        -[b7to1,b7to2,b7to3,b7to4]]/10^9)));
  disp(sprintf('  Less UKI   %7.1f %7.1f %7.1f %7.1f  ', ...
  sumnan([bwto1,bwto2,bwto3,bwto4;
        -[b8to1,b8to2,b8to3,b8to4]]/10^9)));
  disp(sprintf('  Less both  %7.1f %7.1f %7.1f %7.1f  ', ...
  sumnan([bwto1,bwto2,bwto3,bwto4;
        -[b7to1,b7to2,b7to3,b7to4];
        -[b8to1,b8to2,b8to3,b8to4]]/10^9)));
  disp('                                              ')
  disp('Other world totals:                      ')
  disp(sprintf('  IMF        %7.1f %7.1f %7.1f %7.1f  ', ...
  [bito1,bito2,bito3,bito4]/10^9))
  disp(sprintf('  UN         %7.1f %7.1f %7.1f %7.1f  ', ...
  [buto1,buto2,buto3,buto4]/10^9))
  disp(sprintf('  WDI        %7.1f %7.1f %7.1f %7.1f  ', ...
  [bbto1,bbto2,bbto3,bbto4]/10^9))

  data6       = [ [  NaN,b1to2,b1to3,b1to4]/10^9,pop1,gdppc1, ...
                   dfdi1,dfdik1,dgdp1,fdigdp1a*100,fdigdp1b*100;
                  [b2to1,  NaN,b2to3,b2to4]/10^9,pop2,gdppc2, ...
                   dfdi2,dfdik2,dgdp2,fdigdp2a*100,fdigdp2b*100;
                  [b3to1,b3to2,  NaN,b3to4]/10^9,pop3,gdppc3, ...
                   dfdi3,dfdik3,dgdp3,fdigdp3a*100,fdigdp3b*100;
                  [b4to1,b4to2,b4to3,  NaN]/10^9,pop4,gdppc4, ...
                   dfdi4,dfdik4,dgdp4,fdigdp4a*100,fdigdp4b*100;
                  [b5to1,b5to2,b5to3,b5to4]/10^9,pop5,gdppc5, ...
                   dfdi5,dfdik5,dgdp5,fdigdp5*100,NaN;
                  [b6to1,b6to2,b6to3,b6to4]/10^9,pop6,gdppc6, ...
                   dfdi6,dfdik6,dgdp6,fdigdp6*100,NaN];
  Data6(:,j)  = data6(:);

  if icase==1;
    TSto1(j,:)  = [NaN,b2to1,b3to1,b4to1,b5to1,b6to1]/bwto1*100;
    TSto2(j,:)  = [b1to2,NaN,b3to2,b4to2,b5to2,b6to2]/bwto2*100;
    TSto3(j,:)  = [b1to3,b2to3,NaN,b4to3,b5to3,b6to3]/bwto3*100;
    TSto4(j,:)  = [b1to4,b2to4,b3to4,NaN,b5to4,b6to4]/bwto4*100;

    TSfrom4(j,:)= [b1from4,b2from4,b3from4,NaN,b5from4,b6from4]/bwfrom4*100;

  elseif (icase==2)
    TSto1(j,:)  = [NaN,b2to1,b3to1,b4to1,b5to1,b6to1]/ ...
                       sumnan([bwto1;-b7to1;-b8to1])*100;
    TSto2(j,:)  = [b1to2,NaN,b3to2,b4to2,b5to2,b6to2]/ ...
                       sumnan([bwto2;-b7to2;-b8to2])*100;
    TSto3(j,:)  = [b1to3,b2to3,NaN,b4to3,b5to3,b6to3]/ ...
                       sumnan([bwto3;-b7to3;-b8to3])*100;
    TSto4(j,:)  = [b1to4,b2to4,b3to4,NaN,b5to4,b6to4]/ ...
                       sumnan([bwto4;-b7to4;-b8to4])*100;

    TSfrom4(j,:)= [b1from4,b2from4,b3from4,NaN,b5from4,b6from4]/ ...
                       sumnan([bwfrom4;-b7from4;-b8from4])*100;

  else
    TSto1(j,:)  = [NaN,b2to1,b3to1,b4to1,b5to1,b6to1]/ ...
                       sumnan([bwto1;-b7to1])*100;
    TSto2(j,:)  = [b1to2,NaN,b3to2,b4to2,b5to2,b6to2]/ ...
                       sumnan([bwto2;-b7to2;-b8to2])*100;
    TSto3(j,:)  = [b1to3,b2to3,NaN,b4to3,b5to3,b6to3]/bwto3*100;
    TSto4(j,:)  = [b1to4,b2to4,b3to4,NaN,b5to4,b6to4]/ ...
                       sumnan([bwto4;-b8to4])*100;

    TSfrom4(j,:)= [b1from4,b2from4,b3from4,NaN,b5from4,b6from4]/ ...
                       sumnan([bwfrom4;-b8from4])*100;
  end;

  fig1A(j,:)    = [sumnan([b1to4;b2to4;b3to4]),sumnan([bwto4;-b8to4])]/10^9;
  fig1B(j,:)    = [sumnan([b1from4;b2from4;b3from4]),sumnan([bwfrom4;-b8from4])]/10^9;
  fig1C(j,:)    = [sumnan([b2to1;b3to1]),b4to1,bwto1]/10^9;

  TSrpop(j,:)   = [pop1,pop2,pop3,pop4,pop5,pop6];
  TStaup(j,:)   = [taup1,taup2,taup3,taup4,taup5,taup6];
  TSgdppc(j,:)  = [gdppc1,gdppc2,gdppc3,gdppc4,gdppc5,gdppc6];
  TSfdigdp(j,:) = [fdigdp1,fdigdp2,fdigdp3,fdigdp4,fdigdp5,fdigdp6]*100;
  TSpinet(j,:)  = [pinet1,pinet2,pinet3,pinet4,pinet5,pinet6]/10^12;
  TSpingdp(j,:) = [pingdp1,pingdp2,pingdp3,pingdp4,pingdp5,pingdp6]*100;
  TSfditrd(j,:) = [dfdi1,dfdi2,dfdi3,dfdi4,dfdi5,dfdi6]./ ...
                   Data6(49:54,1)'*100;
  TSfdiktrdo(j) = dfdik4o/Data6(52)*100;
  TSfdiktrd(j,:)= [dfdik1,dfdik2,dfdik3,dfdik4,dfdik5,dfdik6]./ ...
                   Data6(49:54,1)'*100;
  TSfdi(j,:)    = [fdigdp1*ngdp1,fdigdp2*ngdp2,fdigdp3*ngdp3, ...
                   fdigdp4*ngdp4,fdigdp5*ngdp5,fdigdp6*ngdp6];
  TSngdp(j,:)   = [ngdp1,ngdp2,ngdp3,ngdp4,ngdp5,ngdp6];
  TStrend(j,:)  = [trend1,trend2,trend3,trend4,trend5,trend6];
  TSdefl(j,:)   = [pgdpd1,pgdpd2,pgdpd3,pgdpd4,pgdpd5,pgdpd6];
end;

disp(' ')
disp(' ')
disp(' ')
disp('=====================================')
disp(sprintf('Inputs for Computer Codes: %4g:%4g',[t0,t1]))
disp('=====================================')
disp(' ')
disp(' ')
disp('    -------------------------------------------------------')
disp('           FDI Inward Flows to Region 1 (as % of World)')
disp('    -------------------------------------------------------')
disp('     From:    EU   Jap-Kor   CHTM     BRI     ROW    Total')
for i=t0:t1
  j           = i-t0+1;
  disp(sprintf('     %4g %7.1f %7.1f %7.1f %7.1f %7.1f  %7.1f', ...
               [i,TSto1(j,2:6),sumnan(TSto1(j,2:6)')])) 
end;
disp(' ')
disp('    -------------------------------------------------------')
disp('           FDI Inward Flows to Region 2 (as % of World)')
disp('    -------------------------------------------------------')
disp('     From:    US   Jap-Kor   CHTM     BRI     ROW    Total')
for i=t0:t1
  j           = i-t0+1;
  disp(sprintf('     %4g %7.1f %7.1f %7.1f %7.1f %7.1f  %7.1f', ...
               [i,TSto2(j,[1,3:6]),sumnan(TSto2(j,[1,3:6])')])) 
end;

disp(' ')
disp('    -------------------------------------------------------')
disp('           FDI Inward Flows to Region 3 (as % of World)')
disp('    -------------------------------------------------------')
disp('     From:    US      EU     CHTM     BRI     ROW    Total')
for i=t0:t1
  j           = i-t0+1;
  disp(sprintf('     %4g %7.1f %7.1f %7.1f %7.1f %7.1f  %7.1f', ...
               [i,TSto3(j,[1:2,4:6]),sumnan(TSto3(j,[1:2,4:6])')])) 
end;

disp(' ')
disp('    -------------------------------------------------------')
disp('           FDI Inward Flows to Region 4 (as % of World)')
disp('    -------------------------------------------------------')
disp('     From:    US      EU   Jap-Kor    BRI     ROW    Total')
for i=t0:t1
  j           = i-t0+1;
  disp(sprintf('     %4g %7.1f %7.1f %7.1f %7.1f %7.1f  %7.1f', ...
               [i,TSto4(j,[1:3,5:6]),sumnan(TSto4(j,[1:3,5:6])')])) 
end;

disp(' ')
disp('    -------------------------------------------------------')
disp('           FDI Outward Flows from Region 4 (as % of World)')
disp('    -------------------------------------------------------')
disp('     From:    US      EU   Jap-Kor    BRI     ROW    Total')
for i=t0:t1
  j           = i-t0+1;
  disp(sprintf('     %4g %7.1f %7.1f %7.1f %7.1f %7.1f  %7.1f', ...
               [i,TSfrom4(j,[1:3,5:6]),sumnan(TSfrom4(j,[1:3,5:6])')])) 
end;

disp(' ')
disp('    -------------------------------------------------------')
disp('           Populations Relative to the United States       ')
disp('    -------------------------------------------------------')
disp('              US      EU   Jap-Kor   CHTM    BRI     ROW  ')
for i=t0:t1
  j           = i-t0+1;
  disp(sprintf('     %4g %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f', ...
               [i,TSrpop(j,:)]))
end;

disp(' ')
disp('    -------------------------------------------------------')
disp('                  Corporate Income Tax Rates               ')
disp('    -------------------------------------------------------')
disp('              US      EU   Jap-Kor   CHTM    BRI     ROW  ')
for i=t0:t1
  j           = i-t0+1;
  disp(sprintf('     %4g %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f', ...
               [i,TStaup(j,:)]))
end;

disp(' ')
disp('    -------------------------------------------------------')
disp('           Per Capita GDP (PPP) Relative to the US         ')
disp('    -------------------------------------------------------')
disp('              US      EU   Jap-Kor   CHTM    BRI     ROW  ')
for i=t0:t1
  j           = i-t0+1;
  disp(sprintf('     %4g %7.1f %7.1f %7.1f %7.1f %7.1f %7.1f', ...
               [i,TSgdppc(j,:)]))
end;

disp(' ')
disp('    -------------------------------------------------------')
disp('                Inward FDI relative to GDP (%)             ')
disp('    -------------------------------------------------------')
disp('              US      EU   Jap-Kor   CHTM    BRI     ROW  ')
for i=t0:t1
  j           = i-t0+1;
  disp(sprintf('     %4g %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f', ...
               [i,TSfdigdp(j,:)]))
end;

disp(' ')
disp('    -------------------------------------------------------')
disp('                   FDI relative to Trend              ')
disp('    -------------------------------------------------------')
disp('              US      EU   Jap-Kor   CHTM    BRI     ROW  ')
for i=t0:t1
  j           = i-t0+1;
  disp(sprintf('     %4g %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f', ...
               [i,TSfditrd(j,:)]))
end;

disp(' ')
disp('    -------------------------------------------------------')
disp('              Cumulated FDI relative to Trend              ')
disp('    -------------------------------------------------------')
disp('              US      EU   Jap-Kor   CHTM    BRI     ROW  ')
for i=t0:t1
  j           = i-t0+1;
  disp(sprintf('     %4g %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f', ...
               [i,TSfdiktrd(j,:)]))
end;

disp(' ')
disp('    -------------------------------------------------------')
disp('             Net portfolio flows relative to GDP (%)    ')
disp('    -------------------------------------------------------')
disp('              US      EU   Jap-Kor   CHTM    BRI     ROW  ')
for i=t0:t1
  j           = i-t0+1;
  disp(sprintf('     %4g %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f', ...
               [i,TSpingdp(j,:)]))
end;

disp(' ')
disp('    -------------------------------------------------------')
disp('             Net portfolio flows     ')
disp('    -------------------------------------------------------')
disp('              US      EU   Jap-Kor   CHTM    BRI     ROW  ')
for i=t0:t1
  j           = i-t0+1;
  disp(sprintf('     %4g %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f', ...
               [i,TSpinet(j,:)]))
end;

figure(1),bar(1990:2010,fig1A(:,2)),hold,bar(1990:2010,fig1A(:,1),'g')
title('Inward FDI to China and Provinces (excluding UKI)')
legend('All','US+EU+JP','Location','Northwest')
axis([1989,2011,0,70])

figure(2),bar(1990:2010,fig1B(:,2)),hold,bar(1990:2010,fig1B(:,1),'g')
title('Outward FDI from China and Provinces (excluding UKI)')
legend('All','US+EU+JP','Location','Northwest')
axis([1989,2011,0,70])

figure(3),bar(1990:2010,fig1C(:,3)),hold, ...
          bar(1990:2010,fig1C(:,1)+fig1C(:,2),'r'), ...
          bar(1990:2010,fig1C(:,1),'g')
title('Inward FDI to United States')
legend('All','CHTM','EU+JP','Location','Northwest')
axis([1989,2011,0,350])

figure(4),plot(1990:2010,TSfdiktrd(:,4),1990:2010,TSfdiktrdo)
title('Cumulative Inward and Outward FDI for China')
legend('Inward','Outward','Location','Northwest')


