clear all
clc

load beliefs.prn

%Tables 4a & 4b:
newdata = [beliefs(:,1:4) zeros(length(beliefs),3) beliefs(:,5:34) zeros(length(beliefs),10)];

newdata(:,5) = newdata(:,3) + (newdata(:,1)-1)*9;
newdata(:,6) = newdata(:,2) - 15 * (newdata(:,4)-1);
newdata(:,7) = newdata(:,6);
p = newdata(find(newdata(:,3) == 1),[1:6 23:37]);

belief = [];

for i = 1:length(newdata)
    a = [newdata(i,:); zeros(15 - newdata(i,6),min(size(newdata)))];
    belief = [belief; a];
end

for i = 1:length(belief)
    if belief(i,1) == 0
        belief(i,[1:6 8:37]) = belief(i-1,[1:6 8:37]);
        belief(i,7) = belief(i-1,7) + 1;
    end
end

for i = 1:length(belief)
    belief(i,38) = belief(i,22 + belief(i,6));
    if belief(i,6) > 1
        belief(i,39) = belief(i,22 + belief(i,6) - 1);
    end
    if belief(i,6) > 2
        belief(i,40) = belief(i,22 + belief(i,6) - 2);
    end
    belief(i,43) = belief(i,7 + belief(i,7));
    if belief(i,7) > belief(i,6)
        belief(i,44) = belief(i,7 + belief(i,7) - 1);
    end
    if belief(i,7) > belief(i,6) + 1
        belief(i,45) = belief(i,7 + belief(i,7) - 2);
    end
    if belief(i,4) > 1
        belief(i,41) = p(find(p(:,1) == belief(i,1) & p(:,6) == belief(i,7) & p(:,4) == belief(i,4) - 1),6 + belief(i,7));
        if belief(i,7) > 1
            belief(i,42) = p(find(p(:,1) == belief(i,1) & p(:,6) == belief(i,7) & p(:,4) == belief(i,4) - 1),5 + belief(i,7));
        end
    end
end

for i = 1:length(belief)
    if belief(i,4) == 1 | belief(i,7) == 1
        belief(i,46) = 0;
    elseif belief(i,6) == belief(i,7)
        belief(i,46) = belief(i,39)*(1 + (belief(i,41) - belief(i,42))/belief(i,42));
    else
        belief(i,46) = belief(i,44)*(1 + (belief(i,41) - belief(i,42))/belief(i,42));
    end
    if belief(i,6) == belief(i,7)
        if belief(i,6) < 3
            belief(i,47) = 0;
        else
            belief(i,47) = belief(i,39)*(1 + (belief(i,39) - belief(i,40))/belief(i,40));
        end
    elseif belief(i,6) == belief(i,7) - 1
        if belief(i,6) < 2
            belief(i,47) = 0;
        else
            belief(i,47) = belief(i,44)*(1 + (belief(i,44) - belief(i,39))/belief(i,39));
        end
    else
        belief(i,47) = belief(i,44)*(1 + (belief(i,44) - belief(i,45))/belief(i,45));
    end
end

for i = 1:length(belief)
    if belief(i,7) < 3 | belief(i,47) == inf
        belief(i,4) = belief(i,4) + 4;
    end
end

belief = [belief dummyvar(belief(:,5))];

y1 = belief(find(belief(:,4) == 1),43); x1 = belief(find(belief(:,4) == 1),47:101); x1(:,37) = [];
y2 = belief(find(belief(:,4) == 2),43); x2 = belief(find(belief(:,4) == 2),46:101); x2(:,38) = [];
y3 = belief(find(belief(:,4) == 3),43); x3 = belief(find(belief(:,4) == 3),46:101); x3(:,38) = [];
y4 = belief(find(belief(:,4) == 4),43); x4 = belief(find(belief(:,4) == 4),46:101); x4(:,[3:11 38]) = [];

[b1,conf1,resid1,bin1,stats1] = regress(y1,x1);
[b2,conf2,resid2,bin2,stats2] = regress(y2,x2);
[b3,conf3,resid3,bin3,stats3] = regress(y3,x3);
[b4,conf4,resid4,bin4,stats4] = regress(y4,x4);

xx1 = inv(x1'*x1);
xx2 = inv(x2'*x2);
xx3 = inv(x3'*x3);
xx4 = inv(x4'*x4);
beta = [[0 b1(1)]; b2(1:2)'; b3(1:2)'; b4(1:2)'];
d(1,:) = [0 sqrt(xx1(1,1)*stats1(4)) stats1(1)];
d(2,:) = [sqrt(xx2(1,1)*stats2(4)) sqrt(xx2(2,2)*stats2(4)) stats2(1)];
d(3,:) = [sqrt(xx3(1,1)*stats3(4)) sqrt(xx3(2,2)*stats3(4)) stats3(1)];
d(4,:) = [sqrt(xx4(1,1)*stats4(4)) sqrt(xx4(2,2)*stats4(4)) stats4(1)];

results1 = [beta d];

'regression results for table 4a'
[results1]

belief(:,47) = 180 - 12*(belief(:,7) - 1);

z1 = belief(find(belief(:,4) == 1),47:101); z1(:,37) = [];
z2 = belief(find(belief(:,4) == 2),47:101); z2(:,37) = [];
z3 = belief(find(belief(:,4) == 3),47:101); z3(:,37) = [];
z4 = belief(find(belief(:,4) == 4),47:101); z4(:,[2:10 37]) = [];

[c1,conf5,resid5,bin5,stats5] = regress(y1,z1);
[c2,conf6,resid6,bin6,stats6] = regress(y2,z2);
[c3,conf7,resid7,bin7,stats7] = regress(y3,z3);
[c4,conf8,resid8,bin8,stats8] = regress(y4,z4);

betac = [c1(1); c2(1); c3(1); c4(1)];

zz1 = inv(z1'*z1);
zz2 = inv(z2'*z2);
zz3 = inv(z3'*z3);
zz4 = inv(z4'*z4);

e(1,:) = [sqrt(zz1(1,1)*stats5(4)) stats5(1)];
e(2,:) = [sqrt(zz2(1,1)*stats6(4)) stats6(1)];
e(3,:) = [sqrt(zz3(1,1)*stats7(4)) stats7(1)];
e(4,:) = [sqrt(zz4(1,1)*stats8(4)) stats8(1)];

results2 = [betac e];

'regression results for table 4b'
[results2]

%Table 5:
round = beliefs(:,2) - 15*(beliefs(:,4) - 1);

for i = 1:length(beliefs)
    blf(i,:) = beliefs(i,4 + round(i));
    price(i,:) = beliefs(i,19 + round(i));
    lagprice(i,:) = beliefs(i,18 + round(i));
end

process = [beliefs(:,1:4) ones(length(beliefs),1) round price-lagprice blf-lagprice];
process = process(find(process(:,6) ~= 1),:);

regset = [];

for i = 1:length(process)
    if process(i,3) == 9 | (process(i,3) == 8 & process(i,1) == 4)
        regset = [regset; [process(i,1:7) mean(process(i+1-process(i,3):i,8))]];
    end
end

for i = 1:4
    x = regset(find(regset(:,4) == i),[5 8]);
    xx = inv(x'*x);
    y = regset(find(regset(:,4) == i),7);
    [r1 r2 r3 r4 r5] = regress(y,x);
    beta(i,:) = r1';
    rho(i) = r5(1);
    stderr(i,:) = sqrt(r5(4)*diag(xx));
end

results3 = [beta stderr rho'];

'regression results for table 5'
[results3]

%Table 6:
pricemat = beliefs(find(beliefs(:,4)>1 & round(:,1) == 15),20:34);
blfmat3 = beliefs(find(beliefs(:,4)>1 & round(:,1) == 3),5:19);
blfmat4 = beliefs(find(beliefs(:,4)>1 & round(:,1) == 4),5:19);

maxblf3 = nanmax(beliefs(find(beliefs(:,4)>1 & round(:,1) == 3),5:19)')';
maxblf4 = nanmax(beliefs(find(beliefs(:,4)>1 & round(:,1) == 4),5:19)')';
maxprice = max(pricemat')';

for i = 1:150
    for j = 1:15
        if pricemat(i,j) == maxprice(i)
            peakprice(i,:) = j;
        end
        if blfmat3(i,j) == maxblf3(i)
            peakblf3(i,:) = j;
        end
        if blfmat4(i,j) == maxblf4(i)
            peakblf4(i,:) = j;
        end
    end
end

dummy2 = beliefs(find(beliefs(:,4)>1 & round(:,1) == 15),4) == 2;
dummy3 = beliefs(find(beliefs(:,4)>1 & round(:,1) == 15),4) == 3;

y1 = peakprice - peakblf3;
y2 = peakprice - peakblf4;
x = [ones(length(dummy2),1) dummy2 dummy3];

[r1,r2,r3,r4,r5] = regress(y1,x);
[q1,q2,q3,q4,q5] = regress(y2,x);
results4 = [3 r1' sqrt(diag(inv(x'*x)*r5(4)))' r5(3);4 q1' sqrt(diag(inv(x'*x)*q5(4)))' q5(3)];

'regression results for table 6'
[results4]

%Tables 7a & 7b:
lagfv = 180 - 12*(round - 2);

for i = 1:length(beliefs)
    sentiment(i,:) = beliefs(i,5 + round(i)) - beliefs(i,4 + round(i));
    price(i,:) = beliefs(i,19 + round(i));
    lagprice(i,:) = beliefs(i,18 + round(i));
    lag2price(i,:) = beliefs(i,17 + round(i));
end

pessimist = sentiment < 0;
process = [beliefs(:,1:4) ones(length(beliefs),1) round price-lagprice lagprice-lagfv lagprice-lag2price];

regset = [];

for i = 1:length(process)
    if process(i,3) == 9 | (process(i,3) == 8 & process(i,1) == 4)
        regset = [regset; [process(i,1:9) sum(pessimist(i+1-process(i,3):i)) mean(sentiment(i+1-process(i,3):i))]];
    end
end
regset = regset(find(regset(:,6) > 2 & regset(:,6) < 15),:);% Eliminating periods 1, 2 and 15.

for i = 1:4
    x = regset(find(regset(:,4) == i),[5 10 8 9]);
    xx = inv(x'*x);
    y = regset(find(regset(:,4) == i),7);
    [r1 r2 r3 r4 r5] = regress(y,x);
    results5(i,:) = [r1' sqrt(r5(4)*diag(xx))' r5(1)];
end
'regression results for table 7a'
[results5]

for i = 1:4
    x = regset(find(regset(:,4) == i),[5 11 8 9]);
    xx = inv(x'*x);
    y = regset(find(regset(:,4) == i),7);
    [r1 r2 r3 r4 r5] = regress(y,x);
    results6(i,:) = [r1' sqrt(r5(4)*diag(xx))' r5(1)];
end
'regression results for table 7b'
[results6]