%URN      Randomly draw probabilities for transition into, out of, and 
%         through the war.  Keep only those that yield time series
%         satisfying the following conditions:
%           1. Probability of being in war state = 15% x range of [.7,1.3]
%           2. Average duration of war = 3.7 years x range of [.7,1.3]
%           3. Probability of entering war = 4.1% x range of [.7,1.3]
%
%         The transition matrix has the form:
%                                  States in t
%                  1          2       3       4       5          6        
%               before     (     states in war            )  ( after war )
%         States
%         in t+1: before   P(1,1)  P(1,2)  ...            P(1,5)   P(1,6)
%                   2      P(2,1)  P(2,2)  ...            P(2,5)   P(2,6)
%                   3      P(3,1)  P(3,2)  ...            P(3,5)   P(3,6)
%                   4      P(4,1)  P(4,2)  ...            P(4,5)   P(4,6)
%                   5      P(5,1)  P(5,2)  ...            P(5,5)   P(5,6)
%                 after    1-sum   1-sum   ...            1-sum    1-sum

%         Ellen McGrattan, 6-22-02
%         Revised, ERM, 1-19-05
%

% User inputs:
ndraw = 500000;   % number of draws
idep  = 2;        % idep = 1 if postwar depression is uncertain
                  %      = 2 if no postwar depression is certain
                  %      = 3 if postwar depression is certain

%___________________________________________________________________________
%disp(sprintf('Drawing %g times from the urn ....',ndraw))
rand('state',1004)
keep  = [];
stats = [];
nkeep = 0;

for m=0:5;
for n=1:ndraw;

  %if (rem(n,1000)==0); 
  %  disp(sprintf('   currently at iteration %6g',n))
  %end;
  P              = zeros(6);
  P(2,1)         = .1*m+.5*rand;
  tem            = rand(5,1);
  P([1,3:6],1)   = tem/sum(tem) * (1-P(2,1));

  P(3,2)         = .1*m+.5*rand;
  tem            = rand(5,1);
  P([1,2,4:6],2) = tem/sum(tem) * (1-P(3,2));

  P(4,3)         = .1*m+.5*rand;
  tem            = rand(5,1);
  P([1:3,5:6],3) = tem/sum(tem) * (1-P(4,3));

  P(5,4)         = .1*m+.5*rand;
  tem            = rand(5,1);
  P([1:4,6],4)   = tem/sum(tem) * (1-P(5,4));

  P(6,5)         = .1*m+.5*rand;
  tem            = rand(5,1);
  P(1:5,5)       = tem/sum(tem) * (1-P(6,5));

  P(6,6)         = .9+.1*rand;
  tem            = rand(5,1);
  P(1:5,6)       = tem/sum(tem) * (1-P(6,6));

  T    = 10000;
  S    = markov(P',T,1);

  i    = find( (S(1:T-11)==1 | S(1:T-11)==6) & (S(2:T-10)>1 & S(2:T-10)<6));

  p1yr = sum(S(i+2)==6)/length(i);
  p2yr = sum(S(i+2)<6 & S(i+3)==6)/length(i);
  p3yr = sum(S(i+2)<6 & S(i+3)<6 & S(i+4)==6)/length(i);
  p4yr = sum(S(i+2)<6 & S(i+3)<6 & S(i+4)<6 & S(i+5)==6)/length(i);
  p5yr = sum(S(i+2)<6 & S(i+3)<6 & S(i+4)<6 & S(i+5)<6 & S(i+5)==6)/length(i);
  p6yr = sum(S(i+2)<6 & S(i+3)<6 & S(i+4)<6 & S(i+5)<6 & S(i+5)<6 & ...
             S(i+6)==6)/length(i);
  p7yr = sum(S(i+2)<6 & S(i+3)<6 & S(i+4)<6 & S(i+5)<6 & S(i+5)<6 & ...
             S(i+6)<6 & S(i+7)==6)/length(i);
  p8yr = sum(S(i+2)<6 & S(i+3)<6 & S(i+4)<6 & S(i+5)<6 & S(i+5)<6 & ...
             S(i+6)<6 & S(i+7)<6 & S(i+8)==6)/length(i);
  p9yr = sum(S(i+2)<6 & S(i+3)<6 & S(i+4)<6 & S(i+5)<6 & S(i+5)<6 & ...
             S(i+6)<6 & S(i+7)<6 & S(i+8)<6 & S(i+9)==6)/length(i);
  p10yr= sum(S(i+2)<6 & S(i+3)<6 & S(i+4)<6 & S(i+5)<6 & S(i+5)<6 & ...
             S(i+6)<6 & S(i+7)<6 & S(i+8)<6 & S(i+9)<6 & S(i+10)==6)/length(i);

  inwar= T-sum(S==6)-sum(S==1);
  pyr  = [p1yr,p2yr,p3yr,p4yr,p5yr,p6yr,p7yr,p8yr,p9yr,p10yr];
  x    = [inwar/T,sum(pyr.*[1:10]),length(i)/T];
  %
  % Check to see if conditions satisfied for keeping draw:
  %
  if ((x(1)>26/171*.7) & (x(1)<26/171*1.3) & (x(2)>3.7*.7) &  ...
      (x(2)<3.7*1.3) & (x(3)>7/171*.7) & (x(3)<7/171*1.3));
    if nkeep>0;
      % 
      % Check for duplication in the list of Markov chains before adding
      % 
      iskip = 0;
      for i=1:nkeep;
        Q = keep((i-1)*6+1:i*6,:);
        s = sum(sum(abs(P-Q)));
        if (abs(s)<.001);
          iskip = 1;
        end;
      end;
      if iskip==0;
        keep   = [keep;P];
        stats  = [stats;x];
        nkeep  = nkeep+1;
        %disp(sprintf('        %4g by iteration %6g',nkeep,n))
      end;
    else
      keep  = P;
      stats = x;
      nkeep = 1;
      %disp(sprintf('  found %4g by iteration %6g',nkeep,n))
    end;     
  end;
end;
end;

Pi = [];
for i=1:nkeep;
  p             = keep((i-1)*6+1:i*6,:);
  pi            = zeros(7);
  pi(1:5,1:5)   = p(1:5,1:5);
  if idep==1;
    pi(6,1)     = .22*p(6,1);
    pi(7,1)     = .78*p(6,1);
    pi(6,2)     = .57*p(6,2);
    pi(7,2)     = .43*p(6,2);
    pi(6,3)     = .53*p(6,3);
    pi(7,3)     = .47*p(6,3);
    pi(6,4)     = .49*p(6,4);
    pi(7,4)     = .51*p(6,4);
    pi(6,5)     = .56*p(6,5);
    pi(7,5)     = .44*p(6,5);
    pi(1:6,6)   = p(:,6);
    pi(7,6)     = 0.;
    pi(1:5,7)   = p(1:5,6);
    pi(7,7)     = p(6,6)-min(.5,p(6,6));
    pi(6,7)     = min(.5,p(6,6));
  end;
  if idep==2;
    pi(1:6,1:6) = p(1:6,1:6);
    pi(7,7)     = 1.;
  end;
  if idep==3;
    pi(7,1:5)   = p(6,1:5);
    pi(1:5,7)   = p(1:5,6);
    pi(6,6)     = 1.;
    pi(7,7)     = p(6,6);
  end;
  Pi   = [Pi;pi];
end;


save urn
%return
diary chains.inp
  disp(nkeep)
  disp(Pi)
diary
