%% Here different agents have different forecasts, different consumption
% and, as a result, different debt or financial wealth.  It is assumed
% that price setting depends only on average forecasts and that income
% (other than interest on debt or savings) is the same for all agents
% There is a limit on allowed debt. If agents hit this limit, they consume
% less reducing aggregate output, inflation, and the interest rate.  



%% Warming up
clear
close all
clc


%% Load shocks
%load matrix_shocks;

%%Generate shocks
%4 
rng('default')
%7 is OK; %9 is OK %12 OK %23 OK
%s = rng;
innovations=randn(5000,1,3);
innovations=randn(5000,1,3);
innovations=randn(5000,1,3);
innovations=randn(5000,1,3);
innovations=randn(5000,1,3);
innovations=randn(5000,1,3);
  innovations=randn(100000,1,3); 
%innovations=randn(5000,1,3);
%innovations=randn(5000,1,3);
% innovations=randn(5000,1,3);
% innovations=randn(5000,1,3);
% innovations=randn(5000,1,3);
% innovations=randn(5000,1,3);
% innovations=randn(5000,1,3);
% innovations=randn(5000,1,3);
 %innovations=randn(5000,1,3);
% innovations=randn(5000,1,3);
% innovations=randn(5000,1,3);
% innovations=randn(5000,1,3);
% innovations=randn(5000,1,3);
% innovations=randn(5000,1,3);
% innovations=randn(5000,1,3);
% innovations=randn(5000,1,3);
% innovations=randn(5000,1,3);
%rng(s);
%innovations=randn(1000,1,3);

%% Parameters of the model
%fundamental parameters
[realiz,rrr,sss]=size(innovations);
bet=0.99;
kappa=0.5;
phi=1;
omega=0.75;
innovations1=innovations(:,1,1);%zeros(realiz,1,1);
innovations2=innovations(:,1,2);
innovations3=innovations(:,1,3);
disturbances=[innovations1,innovations2,innovations3];
mm =10;    %learning parameter gamma BB: intensity of the choice %switching parameter gamma in Brock Hommes
%pstar = 0.02;      % the central bank's inflation target
a1 = 1;%0.5;     %coefficient of expected output in output equation
a2 =-1;%0.20;   %a is the interest elasticity of output demand 
%b1 =0.5;     %b1 is coefficient of expected inflation in inflation equation
b1=(bet/(1+bet*kappa));
% here trouble with monthly vs quarterly b2=0.3 not 0.1
%b2 =0.1275;% b 0.05; %  % is coefficient of output in inflation equation here 0.3 not 0.1 for quarterly not monthly
b2=((1-omega*bet)/(1+bet*kappa))*((1-omega)/omega)*(1+phi)
c1 = 1.5;   %c1 is coefficient of inflation in Taylor equation
c2 = 0.00000001;   %c2 is coefficient of output in Taylor equation
c3 = 0.0;%;%0.5;%;%0.5;%0.7;   %interest smoothing parameter in Taylor equation
b3=(kappa/(1+bet*kappa));      %coef. on lag inflation
A = [1 -b2;-a2*c1 1-a2*c2];
B = [b1 0;-a2 a1]; %Note: p.11 of the book contains a wrong expression for B!
%C = [1-b1 0;0 1-a1]; %Thsi is the  matrix of lagged variables
C = [b3 0;0 1-a1]; %This is the  matrix of lagged variables
smooth = [0;a2*c3]; %This is book's vector b, to multiply lagged interest rate 
 

%T = 10000;
[T,rea]=size(disturbances);

nn=200; % maximum number of iterations in liquidity constrained loop
delta = 0.0000000000000000001; % stop iterating in liquidity constrained loop when change less than delta
sigma1 = 0.0135; %0.0135;   %   %standard deviation shocks output
sigma2 = 0.003;%0.003;%0.005;      %standard deviation shocks inflation
sigma3 = 0.003;%0.005;      %standard deviation shocks Taylor
rho=0.5;           %rho in mean squares errors, BB: memory of agents, see p. 17
rhoout=0.9;           %rho in shocks output
rhoinf=0.8;           %rho in shocks inflation
rhotayl=0.5;      %rho in shocks Taylor
%rhoBH=0.0;          %alpha persistence parameter
%epfs=pstar;           %forecast inflation targeters
rlb =-1; % -.03 lower bound of interest rate - normal interest rate, BB: that's NEW, ok but check:  play with sigma1, sigma2 and sigma3 to make it work
plb =-1; % -.02 lower bound to inflation - target inflation. can be downward nominal rigidity, but really to prevent red dings

maxii=0;
iii = zeros(T,1);
maxzd=0;
%BB: NEW BLOCK
% new individual specific inclinations to extrapolate (identical for
% inflation and output)
N = 100; % number of types of agents 
amax =0.99; %largest fraction extrapolate either inflation or output
amin =0.01; % 1-largest fraction fundamentalists for either inflation or output
dl =1; % debt limit as a share of normal output




%jst here
sever = 0.1;
tlag =4;%120/4;
tlead = 4;%36/4;
ns = 0;
debts = 0;
cdl = 0;
cdld = 0;
ysever=0;




astep = (amax-amin)/(N-2); % fraction of the population in each of the groups which sometimes extrapolates;
pops = ones(1,N)*astep;
pops(1) = amin;
pops(N) = 1-amax; % pops is a row vector with the fraction of the population in each group
debt = zeros(T,1);%for plot
wealth_tot=zeros(T,1);
wealth = zeros(N,1); %distribution of wealth 
wndc = zeros(N,1); %wealth not debt constrained
cdiff = zeros(N,1); 
extyv = zeros(N,1); 
extpv = zeros(N,1);


wextrap = zeros(T,1);
wmiddle = zeros(T,1);
wfund = zeros(T,1);
oops= zeros(T,1);



p = zeros(T,1); %inflation zero for now
y = zeros(T,1); %output gap
plagt = zeros(T,1); %lagged inflation
ylagt = zeros(T,1);%lagged output gap
r = zeros(T,1); %interest rate
rlagt = zeros(T,1);%lagged interest rate 
epf = zeros(T,1);
epc = zeros(T,1);
CRp = zeros(T,1);
FRp = zeros(T,1);
alfapt = zeros(T,1);
eyfunt = zeros(T,1);
CRy = zeros(T,1);
FRy = zeros(T,1);
alfayt = zeros(T,1);
epsilont = zeros(T,1);
etat = zeros(T,1);
ut = zeros(T,1);
zt = zeros(T,1); %effect of liquidity constraint on demand
epsv = zeros(T,1); 
eysv = zeros(T,1);
anspirits=zeros(T,1);
wealth_all=zeros(T,N);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%heuristic model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
    alfap=0.5;
    alfay=0.5;
  
    T=10000;
for t=2:T
%   try
    epsilont(t) = rhoout*epsilont(t-1) + sigma1*disturbances(t,1);    %shocks in output equation (demand shock)
    etat(t)= rhoinf*etat(t-1) + sigma2*disturbances(t,2);      %shocks in inflation equation (supply shock)
    ut(t) = rhotayl*ut(t-1) + sigma3*disturbances(t,3);              %shocks in Taylor rule (interest rate shock)

    
    epsilon = epsilont(t);
    eta = etat(t);
    u = ut(t);

    pstar(t) = 0; %Target inflation 

    epfs = pstar(t);  %Targeters inflation expectation, BB: fundamentalists

    shocks = [eta;a2*u+epsilon]; %shock vector at time t (vector ut in eq. 1.21 in the book)
   
    epcs=p(t-1); %extrapolative inflation expectation
  

    eps=alfap*epcs+(1-alfap)*epfs; %market inflation expectation 
    % if perfect Central Bank credibility set eps=epfs:  all agents are inflation targeters
    epsv(t) = eps; %BB, the vector of market inflation expectation
  
  
    eychar=y(t-1); %extrapolative output gap expectation
    eyfun=0; %fundamentlist output gap expectation (signal is zero here, not in the book though), BB a check?
    eyfunt(t)=eyfun;
    
    eys=alfay*eychar+(1-alfay)*eyfun; %market output gap expectation = alphay*y(t-1)
    %  all agents are fundamentalists set eys=eyfun;
    eysv(t) = eys;  %BB, the vector of output gap expectation
    
    forecast = [eps;eys]; %market inflation and output gap expectation
    
    plag=p(t-1);
    ylag=y(t-1);
    rlag=r(t-1);
    lag = [plag;ylag];
    %vector b in eq. 1.21 in the book
    D = B*forecast + C*lag + smooth*rlag + shocks; %right hand side of eq. 1.21 in the book p. 11
   
    X = A\D;          % inflation and output gap solution (vector) at time t
%         if X==NaN
%         continue
%     end

    y(t) = X(2,1); % output gap if no one is liquidity constrained
    p(t) = X(1,1); % inflation if no one is liquidity constrained
    r(t) = c1*p(t)+c2*y(t)+c3*r(t-1)+u; % interest rate if no one is liquidity constrained
    
    if r(t)<rlb;     
        y(t) = a1*eys+a2*(rlb-eps) + (1-a1)*ylag+epsilon;
        p(t) = b1*eps+b2*y(t) + (1-b1)*plag+eta;
        r(t) = rlb;     
    end;
    
    %p(t) = max(p(t),plb); % inflation can't be less than plb
    
    
    y1 = y(t);
    p1 = p(t);
    r1 = r(t);
    
    
    
    cdiff = a1*(eyfun-eys +extyv*(eychar-eyfun))-a2*(epfs-eps+extpv*(epcs-epfs)); 
    % groups' desired consumption minus average consumption, BB: OK
    
    wealth1= wealth; %last periods ratio of wealth to normal output
    
    %wndc = wealth*exp(r(t-1)-p1) - cdiff; 
     wndc=wealth*(1+r(t-1))*(1-p(t))-cdiff;
    %wndc=wealth*(1+r(t-1))/(1+p(t))-cdiff;
    % groups' desired financial wealth (to avoid crash in the case of inflation lower -100%)
    %definita come wealth disponbile maggiorata degli interessi meno quella
    %necessaria a finanziarie il consumo desiderato
    
    wndc1= wndc;
  %se la ricchezza desiderata è inferiore al vincolo (-1) 
    zv = (wndc<-dl).*(wndc+dl); % vector of effects of liquidity constraint on consumption (differenze tra la mia wealth desiderata e quella effettiva: so it's negative)
    
    
    z = pops*zv; % effect of liquidity constraint on demand interest rate paribus (how much overall consumption non realizzato per via del vincolo)
    for ii=1:nn;
   
    if z==0
        break; % don't waste time adjusting for the effect of the liquidity constraint if it is not binding
    end
 
     Ddc = D+[0;z];%right hand side of eq. 1.21 modified with demand shock due to the liquidity constraint
     zold=z; % save the old effect of the liquidity contraint to stop iterating if it doesn't change
    
     X = A\Ddc; % inflation and output gap solution (vector) at time t
      
     p(t)= X(1,1);  % extracts time t inflation from solution vector
     y(t)= X(2,1); % extracts time t output gap from solution vector
     
     r(t) = c1*p(t)+c2*y(t)+c3*r(t-1)+u;

    if r(t)<rlb;     
        y(t) = a1*eys+a2*(rlb-eps) + (1-a1)*ylag+z+epsilon;
        p(t) = b1*eps+b2*y(t) + (1-b1)*plag+eta;
        r(t) = rlb;     
    end;
       

      
     changey = y(t)-y1; % effect of liquidity constraint on output
     changer = r(t)-r1; % effect of liquidity constraint on the interest rate
    % wndc = wndc1 + changey -a2*changer + wealth1*(exp(r(t-1)))*(exp(-p(t)) - exp(-p1)); % new desired financial wealth given lower income and interest rate
    %wndc = wndc1 + changey -a2*changer + wealth1*(1+r(t-1))*(1/(1+p(t)) - 1/(1+p1));
    wndc = wndc1 + changey - a2*changer + wealth1*(1+r(t-1))*(p1-p(t));
    zv = (wndc<-dl).*(wndc+dl); % vector of effects of liquidity constraint on consumption 
     z = pops*zv; % effect of liquidity constraint on demand updated
%      if z==zold
%          break; % stop iterating if the liquidity constraint effect is the same as last iteration
%      end
      if abs(z-zold)<delta
         maxzd=max(maxzd,abs(z-zold));
         break; % stop iterating if the liquidity constraint effect is the same as last iteration
     end

    end
 maxii = max(ii,maxii);
    iii(t) = ii; %store how many iterations were needed to get to a liquidity constraint fixed point
    zt(t) = z; %store liquidity constraint effect
    
    wealth = wndc-zv; % wealth with borrowing limit 
    
    
    
    
    debtv = -wealth.*(wealth<0);
    debt(t) = pops*debtv;
        % jst here
    if t>tlag+tlead
    if y(t-tlead)<-sever & y(t-tlead-1)>-sever
        ns = ns+1;
        debts = debts+debt(t-tlead);
        cdl = cdl + debt(t-tlead)-debt(t-tlag-tlead);
        cdld = cdld + debt(t)-debt(t-tlead);
      
    end
    end   
    wealth_totv=wealth.*(wealth>0);
     wealth_tot(t)=pops* wealth_totv;
    wextrap(t) = wealth(1); %  financial wealth of those who always extrapolate
    wmiddle(t) = wealth(floor(N/2)); % financial wealth of a middling group
    wfund(t) = wealth(N); % financial wealth of those who never extrapolate
    oops(t) = pops*wealth; % total financial wealth which should be zero
        wealth_all(t,:)=wealth;
    
    
    %So far for time t variables, now we need to construct the alphas and
    %betas for next period: individuals revise their decision rules based on the
    %discrepancies between period t variables and their previous period
    %expectations of them
    
    plagt(t)=p(t-1);
    ylagt(t)=y(t-1);   
    rlagt(t)=r(t-1);   
    
    
    CRp(t) = rho*CRp(t-1) - (1-rho)*(epcs-p(t))^2;
    FRp(t) = rho*FRp(t-1) - (1-rho)*(epfs-p(t))^2;
    CRy(t) = rho*CRy(t-1) - (1-rho)*(eychar-y(t))^2;
    FRy(t) = rho*FRy(t-1) - (1-rho)*(eyfun-y(t))^2;
  

    
    
    % here extreme y causes Matlab to get to 0/0 so to prevent that 
    badp  = (exp(mm*CRp(t))+exp(mm*FRp(t)))==0;
    bady = (exp(mm * CRy(t)) + exp(mm*FRy(t)))==0;
    
    
    alfap =(amax*badp+exp(mm*CRp(t)))/(badp+exp(mm * CRp(t)) + exp(mm*FRp(t))); %more generally, with persistence in alpha: rhoBH*alfapt(t-1)+(1-rhoBH)*exp(mm*CRp(t))/(exp(mm * CRp(t)) + exp(mm * FRp(t)));
    alfay = (amax*bady+exp(mm*CRy(t)))/(bady+exp(mm * CRy(t)) + exp(mm*FRy(t))); %more generally, with persistence in alpha: rhoBH*alfayt(t-1)+(1-rhoBH)*exp(mm*CRy(t))/(exp(mm * CRy(t)) + exp(mm * FRy(t)));
   
   % almost all are inflation fundamentalists (extremely credible inflation target)
   % alfap = amin;
   % almost all are output fundamentalists
   % alfay = amin;
    
    nextp = 1+(floor((alfap-amin)/astep))*(alfap>amin); % the number of groups  who extrapolate inflation
    nextp = nextp + (nextp>(N-1))*(N-1-nextp); % alphap may have to be be rounded down to the allowed maximum alfam  
    nexty = 1+floor((alfay-amin)/astep)*(alfay>amin); % the number of groups  who extrapolate output
    nexty = nexty + (nexty>(N-1))*(N-1-nexty); % alphap may have to be rounded down to the allowed maximum alfam
    
    extpv=zeros(N,1);
    extpv(1:(nextp)) = 1; % which groups extrapolate inflation
  

   
    extyv = zeros(N,1);
    extyv(1:(nexty)) = 1; % which groups extrapolate output
    
    alfap = amin+astep*(nextp-1); % round alfap to the minimum alfam plus an even number of allowed steps.a finite number of types
    alfay = amin+astep*(nexty-1);
    
     
    
    alfapt(t) = alfap;
    alfayt(t) = alfay;

%animal spirits see p. 12
    if eychar>0;
        anspirits(t)=alfay;
    end
    if eychar<0;
        anspirits(t)=1-alfay;
    end
   

end

%jst here
debts = debts/ns;
cdl = cdl/ns;
cdld = cdld/ns;
'mean, median, max, min, standard deviation, kurtosis, & skewness of the output gap'


Mean_p    = mean(p)
Stadev_p  = std(p)
Median_p  = median(p)
Max_p     = max(p)
Min_p     = min(p)


Mean_r    = mean(r)
Stadev_r  = std(r)
Median_r  = median(r)
Max_r     = max(r)
Min_r     = min(r)



% %%%%%%%%%%%%%%%%%%%%%%%%%
tstart= 0.1*T; %drop the first 10% observations
%T=4900;
dy = (y(tstart:T)-y(tstart-1:T-1))*100;
dp=(p(tstart:T)-p(tstart-1:T-1))*100;
rr = (r(tstart-1:T-1)-(p(tstart:T)))*100;
yl=(y(tstart-1:T-1))*100;
pl = (p(tstart-1:T-1))*100;
rl = (r(tstart-1:T-1))*100;
y = (y(tstart:T))*100;
p = (p(tstart:T))*100;
r = (r(tstart:T))*100;
anspiritsl=(anspirits(tstart:T))*100;
epsvl = (epsv(tstart-1:T-1))*100;
eysvl = (eysv(tstart-1:T-1))*100;


alfayt = alfayt(tstart:T);
alfapt = alfapt(tstart:T);

wfund = (wfund(tstart:T))*100;  % expressed as % of Y
wextrap = (wextrap(tstart:T))*100; % % expressed as % of Y
wmiddle=(wmiddle(tstart:T))*100;
zt = (zt(tstart:T))*100; %effects of LC on aggregate demand
debt =(debt(tstart:T))*100;  % expressed as % of Y
wealth_tot =(wealth_tot(tstart:T))*100;  % expressed as % of Y
wealth_all =(wealth_all(tstart:T,:))*100;




Mean_p    = mean(p)
Stadev_p  = std(p)
Median_p  = median(p)
Max_p     = max(p)
Min_p     = min(p)



Mean_r    = mean(r)
Stadev_r  = std(r);
Median_r  = median(r)
Max_r     = max(r)
Min_r     = min(r)


Mean_y    = mean(y)
Stadev_y  = std(y)
Median_y  = median(y)
Max_y    = max(y)
Min_y    = min(y)




 %jst here
 ' if output less than sever then debt, change in debt over past tlag periods change in debt over next tlead periods'
 [debts,cdl,cdld]

figure('position', [100, 100, 700,250]); clf; hold on;
figure(1); clf; hold on;
subplot(1,2,1)
scatter(debt,y,'.');
grid on;
xlabel('Debt','interpreter', 'latex');
ylabel('Output','interpreter', 'latex')
title({'Output and debt'},'interpreter', 'latex','FontSize',11)
subplot(1,2,2)
scatter(zt,y,'.');
grid on;
xlabel('Liquidity constraint effect','interpreter', 'latex');
ylabel('Output','interpreter', 'latex');
title({'Output and debt limit effect'}, 'interpreter', 'latex','FontSize',11);


figure(13); clf; hold on;
plot(tstart:T,wextrap);
grid on;
xlabel('Time');
ylabel('wealth of extrapolators');
title({'Wealth of Those who Always extrapolate'},'FontWeight','bold');



figure(14); clf; hold on;
plot(tstart:T,wfund);
grid on;
xlabel('Time');
ylabel('wealth of fundamentalists');
title({'Wealth of Those who Never Extrapolate'},'FontWeight','bold');


%figure(15); clf; hold on;
%plot(1:T,oops);
%grid on;
%xlabel('Time');
%ylabel('total wealth');
%title({'total wealth'},'FontWeight','bold');

figure(18); clf; hold on;
plot(tstart:T,zt);
grid on;
xlabel('Time');
ylabel('liquidity constraint effect')
title({'Reduction in demand due to liquidity constraint'},'FontWeight','bold');

figure(19); clf; hold on;
plot(zt,y);
grid on;
xlabel('liquidity constraint effect');
ylabel('Output Gap');
title({'Output gap and liquidity constraint'},'FontWeight','bold');


figure(21); clf; hold on;
plot(debt,y);
grid on;
xlabel('debt');
ylabel('Output Gap')
title({'Output and Debt'},'FontWeight','bold')

figure(137); clf; hold on;
plot(r,y);
grid on;
xlabel('r');
ylabel('y');
title({'Interest rate & output'},'FontWeight','bold');

%figure(16); clf; hold on;
%plot(wfund,y);
%grid on;
%xlabel('wealth of fundamentalists');
%ylabel('output gap');
%title({'output and wealth of fundamentalists'},'FontWeight','bold');

%figure(17); clf; hold on;
%plot(wextrap,y);
%grid on;
%xlabel('wealth of extrapolators');
%ylabel('output gap');
%title({'output and wealth of extrapolators'},'FontWeight','bold');

figure(1);
plot(tstart:T,y);
grid on;
xlabel('Time');
ylabel('Level');
title({'output'},'FontWeight','bold');