%% Warming up
clear
close all
clc

%% Parameters of the model


mm = 1;    %learning parameter gamma
%pstar = 0.02;      % the central bank's inflation target
a1 = 1;     %coefficient of expected output in output equation
a2 = -0.2;   %a is the interest elasticity of output demand
b1 = 1;     %b1 is coefficient of expected inflation in inflation equation
b2 = 0.05;   %b2 is coefficient of output in inflation equation
c1 = 1.5;   %c1 is coefficient of inflation in Taylor equation
c2 = 0.5;   %c2 is coefficient of output in Taylor equation
c3 = 0.0;   %interest smoothing parameter in Taylor equation
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 book's matrix of lagged variables
smooth = [0;a2*c3]; %This is book's vector b, to multiply lagged interest rate 
 
T = 12000;
sigma1 = 0.05;      %standard deviation shocks output
sigma2 = 0.05;      %standard deviation shocks inflation
sigma3 = 0.05;       %standard deviation shocks Taylor
rho=0.5;           %rho in mean squares errors
rhoout=0.0;           %rho in shocks output
rhoinf=0.0;           %rho in shocks inflation
rhotayl=0.0;      %rho in shocks Taylor
rhoBH=0.0;          %alpha persistence parameter
%epfs=pstar;           %forecast inflation targeters

p = zeros(T,1); %inflation
pstar = zeros(T,1); %inflation
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
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);
anspirits = zeros(T,1);
epsilont = zeros(T,1);
etat = zeros(T,1);
ut = zeros(T,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%heuristic model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
    alfap=0.5;
    alfay=0.5;
  
for t=2:T
    
    epsilont(t) = rhoout*epsilont(t-1) + sigma1*randn;    %shocks in output equation (demand shock)
    etat(t)= rhoinf*etat(t-1) + sigma2*randn;             %shocks in inflation equation (supply shock)
    ut(t) = rhotayl*ut(t-1) + sigma3*randn;               %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

    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
  
  
    eychar=y(t-1); %extrapolative output gap expectation
    eyfun=0; %fundamentlist output gap expectation (signal is zero here, not in the book though)
    eyfunt(t)=eyfun;
    
   eys=alfay*eychar+(1-alfay)*eyfun; %market output gap expectation = alphay*y(t-1)
   
     %  all agents are fundamentalists set eys=eyfun;
    
    
    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
    X = A\D;          % 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; %obtains time t interest rate from Taylor rule
    
    %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);   
    
    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;
  

    alfap =exp(mm*CRp(t))/(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 = exp(mm*CRy(t))/(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)));

    alfapt(t) = alfap;
    alfayt(t) = alfay;
    
    if eychar>0;
        anspirits(t)=alfay; %this is when extrapolators forecast positive output gap
    end
    if eychar<0;
        anspirits(t)=1-alfay; %this is when extrapolators forecast negative output gap: in this case the "optimists" are the fundamentalists
    end
    

end

autocory = corrcoef(y,ylagt)
autocorp = corrcoef(p,plagt)
coroutputanimal = corr(y,anspirits)
coroutputepsilon = corr(y,epsilont)
coroutputeta = corr(y,etat)
coroutputu = corr(y,ut)
corInflanimal = corr(p,anspirits)
corInflepsilon = corr(p,epsilont)
corInfleta = corr(p,etat)
corInflu = corr(p,ut)

%%  mean, median, max, min, standard deviation, kurtosis

Mean    = mean(y);
Median  = median(y);
Max     = max(y);
Min     = min(y);
Stadev  = std(y);
Kurt    = kurtosis(y);
Skew    = skewness(y);

%% jarque-bera test

[jb,pvalue,jbstat] = jbtest(y,0.05) %test for normality of distribution of output gaps

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tstart= 1001; %drop the first 1000 observations
y = y(tstart:T);
p = p(tstart:T);
r = r(tstart:T);

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


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

figure(100);
hist(y,50); %bins the elements of Y into 50 equally spaced containers
title({'histogram output gap'},'FontWeight','bold');

figure(2); clf; hold on;
plot(tstart:T,p);
grid on;
xlabel('Time');
ylabel('Level');
title({'inflation'},'FontWeight','bold');


figure(200);
hist(p,50); %bins the elements of p into 50 equally spaced containers
title({'histogram inflation'},'FontWeight','bold');

figure(3); clf; hold on;
plot(tstart:T,r);
grid on;
xlabel('Time');
ylabel('Level');
title({'interest rate'},'FontWeight','bold');

figure(201);
hist(r,50); %bins the elements of p into 50 equally spaced containers
title({'histogram interest rate'},'FontWeight','bold');

figure(40); clf; hold on;
plot(tstart:T,anspirits);
grid on;
xlabel('Time');
ylabel('Level');
title({'animal spirits'},'FontWeight','bold');

figure(400);
hist(anspirits,50);
title({'histogram animal spirits'},'FontWeight','bold');

figure(4); clf; hold on;
plot(tstart:T,alfayt);
grid on;
xlabel('Time');
ylabel('Level');
title({'fraction output extrapolators'},'FontWeight','bold');

figure(5); clf; hold on;
plot(tstart:T,alfapt);
grid on;
xlabel('Time');
ylabel('Level');
title({'fraction inflation extrapolators'},'FontWeight','bold');

figure(402);
hist(alfapt,50);
title({'histogram inflation extrapolators'},'FontWeight','bold');

figure(6); clf; hold on;
subplot(4,1,1),plot(tstart:T,y), title('Output gap');
subplot(4,1,2),plot(tstart:T,anspirits), title('Animal Spirits');
subplot(4,1,3),plot(tstart:T,p), title('Inflation');
subplot(4,1,4),plot(tstart:T,r), title('Interest rate');

figure(9); clf; hold on;
subplot(4,1,1),plot(tstart:T,p), title('Inflation');
subplot(4,1,2),plot(tstart:T,alfapt), title('Fraction Inflation Extrapolators');
subplot(4,1,3),plot(tstart:T,r), title('Interest rate');
subplot(4,1,4),plot(tstart:T,y), title('Output Gap');

figure(7); clf; hold on;
plot(anspirits,y);
grid on;
xlabel('Animal Spirits');
ylabel('Output');
title({'Animal spirits and Output Gap'},'FontWeight','bold');

figure(8); clf; hold on;
plot(anspirits,p);
grid on;
xlabel('Animal Spirits');
ylabel('Inflation');
title({'Animal spirits and Inflation'},'FontWeight','bold');

figure(12); clf; hold on;
plot(alfapt,y);
grid on;
xlabel('Extrapolators');
ylabel('Output');
title({'Extrapolators and Output Gap'},'FontWeight','bold');

figure(10); clf; hold on;
plot(alfapt,p);
grid on;
xlabel('Extrapolators');
ylabel('Inflation');
title({'Extrapolators and Inflation'},'FontWeight','bold');

figure(11); clf; hold on;
plot(alfapt,r);
grid on;
xlabel('Extrapolators');
ylabel('Interest rate');
title({'Extrapolators and Interest rate'},'FontWeight','bold');