%Analysis of a stylized technology-only rbc model, here I implement
%a transformation and Metropolis-Hastings method as well
close all;  clear all;  clc;
addpath(genpath('C:/Matlab7/toolbox/KalmanAll'))
%data
[time, gdp, cons, inv1, hours,inv2] = textread('datarbc.txt','%d %f %f %f %f %f');
T=size(cons,1);
gdphp=hpfilter(gdp,100);    y0=log(gdp./gdphp);    
conshp=hpfilter(cons,100);  c0=log(cons./conshp);
hourshp=hpfilter(hours,100);    l0=hours./hourshp-1;%l0=log(hours./hourshp);

[capital] = textread('capital.txt','%f');
capitalhp=hpfilter(capital,100);    k0=log(capital./capitalhp);    

%predefined parameters
beta=0.96;      g=0.0018;        ro=3;
alpha=1/3;      fi=0.17;         lambda=0.65;    delta=0.1;
sigmay=0.05;     sigmak=0.0000000001;     sigmac=0.0001;    sigmal=0.00015;
vol=0.00007;
%set of parameters
q0=[beta; g; ro; alpha; fi; lambda; delta; sigmay; sigmak; sigmac; sigmal; vol];
qmax=[1; 1; 10; 1; 1; 1; 1; 1; 1; 1; 1; 1];
qmin=[0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0];
q1=transform(q0,qmax,qmin,1);
y2=[y0,k0,c0,l0];

% find likelihood maximizer with csminwel - Classical inference
%H0=0.0001*eye(12);  crit=0.00001;   nit=1000;
%[fh,xh,gh,H,itct,fcount,retcodeh] = csminwel(@neolike2, q1, H0, [], crit, nit, qmax, qmin, y2);

%use Metropolis-Hastings + Bayesian inference
%0) set prior variances:
a=0.1; st=0.5; vo=1;
s2=[a; st; a; a; st; st; st; vo; vo; vo; vo; vo];

%set initial values and compute initial likelihood and probability
M=30000;
z(:,1)=q1;
L0=exp(neolike2(z(:,1),qmax, qmin, y2));
P0=normpdf(z(:,1),q1,s2);
P0=exp(log(P0)'*ones(size(P0)));
rej=0;
%1) generate 
for i=2:M
    z(:,i)=z(:,i-1)+s2.*randn(size(s2))/460;
    %new likelihood
    L1=exp(neolike2(z(:,i),qmax, qmin, y2));
    P1=normpdf(z(:,i),q1,s2);
    P1=exp(log(P1)'*ones(size(P1)));
    ratio=L1*P1/(L0*P0);
    anton=rand;
    % accept/reject rule
    if anton<ratio
        L0=L1;  P0=P1;
    else
        z(:,i)=z(:,i-1);
        rej=rej+1;
    end
end

%rejection rate - should be between 25% and 40%
rej/M

%omit the first 1000 observations as 
z=z(:,1000:size(z,2));
%z=z(:,int16(size(z,2)/3):size(z,2));
M=size(z,2);

%transform the results back to orignal parameter values
z1=transform(z, qmax*ones(1,M), qmin*ones(1,M),2);

%plot the histograms of empirical distributions of the parameters
figure(1)
for i=1:12
    subplot(3,4,i)
    xin=z1(i,:); nbins=50;
    xin = reshape( xin, numel(xin), 1 );
    minXin = min(xin); maxXin = max(xin);
    if minXin == maxXin
         bar(minXin,1);
        axis([minXin - 10, minXin + 10, 0, 1]);
    else
         step = (maxXin - minXin) / (nbins-1);
        binc = minXin : step : maxXin;     
        [N, X] = hist(xin, binc);
        bar(X, N/sum(N));
    end
    if i==1
        title('\beta');
    elseif i==2
        title('\gamma');
    elseif i==3
        title('\rho');
    elseif i==4
        title('\alpha');
    elseif i==5
        title('\phi');
    elseif i==6
        title('\lambda');
    elseif i==7
        title('\delta');
    elseif i==8
        title('\sigma_{y}');
    elseif i==9
        title('\sigma_{k}');
    elseif i==10
        title('\sigma_{c}');
    elseif i==11
        title('\sigma_{l}');
    elseif i==12
        title('\sigma_{a}');
    end
    grid on;
end

%compute the resulting means and variances for posterior (zm) and prior (z0)
zm=median(z1')';
zs=sqrt(mean(((z1-zm*ones(1,M)).^2)')');
zm0=q0;
zs0=q0-transform(q1+s2,qmax,qmin,2);

%plot the stylized graphs of priors and posteriors
figure(2)
for i=1:12
    subplot(3,4,i)
    t=max(zm(i)-2.5*zs(i),qmin(i)):(zs(i)/100):min(zm(i)+2.5*zs(i),qmax(i));
    t_=max(zm(i)-5*zs(i),qmin(i)):(zs(i)/100):min(zm(i)+5*zs(i),qmax(i));
    d1=normpdf(t,zm(i),zs(i));
    t2=(transform(t_,qmax(i),qmin(i),1));
    d2=normpdf(t2,q1(i),s2(i));
    d1=d1/max(d1);
    d2=d2/max(d2);
    plot(t,d1,'b',t_,d2,'-r',...
        'LineWidth',2);
    if i==5
        legend('posterior','prior');
    end
    if i==1
        title('\beta');
    elseif i==2
        title('\gamma');
    elseif i==3
        title('\rho');
    elseif i==4
        title('\alpha');
    elseif i==5
        title('\phi');
    elseif i==6
        title('\lambda');
    elseif i==7
        title('\delta');
    elseif i==8
        title('\sigma_{y}');
    elseif i==9
        title('\sigma_{k}');
    elseif i==10
        title('\sigma_{c}');
    elseif i==11
        title('\sigma_{l}');
    elseif i==12
        title('\sigma_{a}');
    end
end

%name the median posteriors
xh=zm;
beta=xh(1); g=xh(2);    ro=xh(3);   alpha=xh(4);    fi=xh(5);   lambda=xh(6);
delta=xh(7);    sigmay=xh(8);   sigmak=xh(9);   sigmac=xh(10);  sigmal=xh(11);  vol=xh(12);
[beta, g, ro, alpha, fi ,lambda, delta, sigmay, sigmak, sigmac, sigmal, vol]
%%%%%%%%%%%%%%%%%%%
%resolve the model
%new parameters - definitions
teta=((((1+g)^(1-fi*(1-ro)))/beta)-(1-delta))/alpha;%%=cap/gdp
nu=1-(delta+g)/teta; %=cons/gdp
labor=1-1/(1+fi*(1-alpha)/(nu*(1-fi)));
%steady state values
cap=labor/(teta^(1/(1-alpha)));
gdp=teta*cap;
cons=nu*gdp;
%more parameter definitions
phi=alpha*beta*gdp/((1+g)*cap);
omega=fi*(1-ro)-1;
a1=labor/(1-labor);
pi=(1-fi)*(1-ro)*a1;
%matrices
D=[ 0,-phi,-omega,pi,0,phi,omega,-pi;
    -gdp,(1+g)*cap,cons,0,0,0,0,0;
    1,0,0,-(1-alpha),-1,0,0,0;
    1,0,-1,-(1+a1),0,0,0,0;
    0,0,0,0,1,0,0,0;
    1,0,0,0,0,0,0,0;
    0,1,0,0,0,0,0,0;
    0,0,1,0,0,0,0,0;];
E=[ 0,0,0,0,0,0,0,0;
    0,(1-delta)*cap,0,0,0,0,0,0;
    0,alpha,0,0,0,0,0,0;
    0,0,0,0,0,0,0,0;
    0,0,0,0,lambda,0,0,0;
    0,0,0,0,0,1,0,0;
    0,0,0,0,0,0,1,0;
    0,0,0,0,0,0,0,1;];
%use gensys
[G,C,impact,fmat,fwt,ywt,gev,eu]=gensys(D,E,zeros(8,1),[0;0;0;0;1;0;0;0],[0,0,0; 0,0,0; 0,0,0; 0,0,0; 0,0,0; 1,0,0; 0,1,0; 0,0,1]);

F = G([1,2,3,4,5],[1,2,3,4,5]);
%H = [1 0 0 0 0; 0 0 1 0 0; 0 0 0 1 0];
H = [eye(4),zeros(4,1)];%G([1,2,3,4],[2,5]);
Q = vol*impact([1,2,3,4,5])*impact([1,2,3,4,5])';
%R = diag([sigmay,sigmac,sigmal]);
R = diag([sigmay,sigmak,sigmac,sigmal]);
initx = [y0(1);k0(1);c0(1);l0(1);0];
%initx = [y0(1);0;c0(1);l0(1);0];
initV = 0.001*eye(5);

%data for the kalman filter
%y1=[y0,c0,l0]';    
y1=y2';
%filter the data
[xfilt, Vfilt, VVfilt, loglik] = kalman_filter(y1, F, H, Q, R, initx, initV);
%smoothing is not needed

%plot the filtered date
y1=y2';
figure(3)
clf
for i=1:4
    subplot(2,2,i)
    hold on
    plot(time, y1(i,:), 'b*-');
    errorbar(time, xfilt(i,:), sqrt(reshape(Vfilt(i, i, :),1,size(Vfilt,3))), 'rx:');
    hold off
    legend('observed', 'filtered', 2)
    xlabel('time')
    if i==1
        ylabel('output')
    elseif i==2
        ylabel('capital')
    elseif i==3
        ylabel('consumption')
    elseif i==4
        ylabel('labor')
    end
end

%compare output with the recovered "tfp shock"
figure(4)
clf
hold on
plot(time, y1(1,:), 'b*-');
errorbar(time, xfilt(5,:), sqrt(reshape(Vfilt(5, 5, :),1,size(Vfilt,3))), 'rx:');
hold off
legend('output', 'productivity filtered', 2)
xlabel('time')
ylabel('productivity vs output')
