%COMPSTEADY    compute the steady states for the McGrattan-Waddle
%              BREXIT model.
%

%              Ellen McGrattan, 4-26-18

%
% Initial Guesses for A,Sig,Ice and steady state for I=5 model 
% with countries:
%
%    1) US+Canada
%    2) EU
%    3) UK 
%    4) Norway+Switzerland
%    5) Japan+Korea+China 
%
% and J=0,1,2 options on trade:
%
%    0) No costs
%    1) Set costs to match goods trade
%    2) Set costs to match goods and services trade
%

I = 5;
J = 1; 

if (J==0); icase=0; end;
if (J==1); icase=1; end;
if (J==2); icase=2; end;

i = int2str(icase);
j = int2str(J);

eval(['load ./cases/xguess',i,'.dat;']);
eval(['load ./cases/iparam',i,'.dat;']);
eval(['load ./cases/moments',j,'.dat;']);
eval(['xguess  = xguess',i,';']);
eval(['iparam  = iparam',i,';']);
eval(['moments = moments',j,';']);

%
% Put moments of interest in datastats
%

datastats    = moments;

%
% If problems with convergence, start with model prediction:
%

%[r,modstats]= steady(xguess,iparam);
%datastats   = modstats;


%
% Indices in iparam for A,sig,ice:
%
iparind      = [14:13+I, ...
                14+11*I+3*I*I:13+11*I+4*I*I, ...
                14+11*I+2*I*I:13+11*I+3*I*I]';

%
% Impose that UK in EU
%
ii           = iparind([13,17,38,42]);
iparam(ii)   = [1;1;0;0];

%
% Indices in data statistics to change:
%

iA           = [2:I]';              % indicies in iparind for A(2:I)
isig1        = [7:10]';             % indicies in iparind for sig(:,1)
isig2        = [11,13:15]';         % indicies in iparind for sig(:,2)
isig3        = [16:17,19:20]';      % indicies in iparind for sig(:,3)
isig4        = [21:23,25]';         % indicies in iparind for sig(:,4)
isig5        = [26:29]';            % indicies in iparind for sig(:,5)
iice1        = [32:35]';            % indicies in iparind for ice(:,1)
iice2        = [36,38:40]';         % indicies in iparind for ice(:,2) 
iice3        = [41:42,44:45]';      % indicies in iparind for ice(:,3)
iice4        = [46:48,50]';         % indicies in iparind for ice(:,4)
iice5        = [51:54]';            % indicies in iparind for ice(:,5)

if (J==0);
  ii         = iparind(31:55);
  iparam(ii) = 0;
  dataind    = [iA;isig1;isig2([1,3,4]);isig3([1,3,4]);isig4;isig5];
end;
if (J==1);
  dataind    = [iA;isig1;isig2([1,3,4]);isig3([1,3,4]);isig4;isig5;
               [32,33,36,41,45,46:48,50:54]'];
end;
if (J==2);
  dataind    = [iA;isig1;isig2([1,3,4]);isig3([1,3,4]);isig4;isig5;
               [33,36,40,41,44,45,46:48,50:54]'];
end;

nind         = length(dataind);


%
% Find new parameters and associated steady state:
% (Some tips: 
%    (1) compute steady state (in steady.m) for initial guess of A,sig
%    (2) use model statistics (modstats) in place of datastats initially
%    (3) homotopy between these modstats and actual datastats
%    (4) subsequently add ice using the additional lines commented below)
%

fullpar      = [I;xguess;iparam;datastats;iparind;nind;dataind];
oldpar       = iparam(iparind(dataind));
newpar       = secant('callsteady',oldpar,fullpar,1e-5,100);

[r,xnew]     = callsteady(newpar,fullpar);
pnew         = iparam;
ii           = iparind(dataind);
pnew(ii)     = newpar;

[r,modstats] = steady(xnew,pnew);

xguess       = xnew;
iparam       = pnew;
tab          = [[1:55]',iparam(iparind),datastats,modstats,datastats-modstats];

disp('-------------------------------------------------')
disp('              Steady state results')
disp('-------------------------------------------------')
disp('   #   Parameter  Moment     Model    Difference ')
disp('-------------------------------------------------')
for i=1:55;
  disp(sprintf(' %3g    %7.4f  %8.4f   %8.4f   %8.4f',tab(i,:)));
end;


return

%
% If iceberg costs hit corners, gradually update costs by running a loop over:
%
iceind       = [iice1;iice2;iice3;iice4;iice5];
mice         = modstats(iceind);
dice         = datastats(iceind);
perdev       = abs((mice-dice)./dice);
iparam(iparind(iceind)) = iparam(iparind(iceind)).* ...
         (1+.01*sign(mice-dice).*perdev);

