% =======================================================================
% This script estimates a simplification of the UCSV-MA model used in
%
% Chan, J.C.C. (2013). Moving Average Stochastic Volatility Models 
% with Application to Inflation Forecast, Journal of Econometrics, 
% 176 (2), 162-172.
%
% The simplifications used are:
%   the MA-coefficient psi is set to 0
%   the AR-coefficient phih is set to 1
%   the AR-coefficient phig is set to 1
%   the constant muh is set to 0
%   the constant mug is set to 0
% such that the only parameters to be estimated are sigh2 and sigg2
% 
% y(t)     = tau(t) + eps_y(t);
% eps_y(t) = u(t) + 0 * u(t-1);                       u(t) ~ N(0,exp(h(t))
% tau(t) = tau(t-1) + eps_tau(t);                       eps_tau(t) ~ N(0,exp(g(t))
% h(t) = 0 + 1*( h(t-1) - 0 ) + eps_h(t)      eps_h(t) ~ N(0,sig_h^2)
% g(t) = 0 + 1*( g(t-1) - 0 ) + eps_g(t)      eps_g(t) ~ N(0,sig_g^2)
%
% =======================================================================

%% prior
%invVpsi = 1;
phih0 = 1; invVphih = 999;
phig0 = 1; invVphig = 999;
muh0 = 0; invVmuh = 999;
mug0 = 0; invVmug = 999;
nuh0 = 10; Sh0 = .05*(nuh0-1);
nug0 = 10; Sg0 = .05*(nug0-1);

disp('Starting UCSV-noMA.... ');
disp(' ' );
start_time = clock;    
    
% initialize the Markov chain
sigh2 = .05; sigg2 = .05;
phih = 1;  phig = 1;
muh = 0; mug = 0;
h = log(var(y)*.5)*ones(T,1);
g = log(var(y)*.5)*ones(T,1);
%psihat = 0;
%psi = psihat;
H = speye(T) - sparse(2:T,1:(T-1),ones(1,T-1),T,T);
%Hpsi = speye(T) + sparse(2:T,1:(T-1),psi*ones(1,T-1),T,T); 
 Hpsi = speye(T) + sparse(2:T,1:(T-1), 0 *ones(1,T-1),T,T); 
%invDpsic = .01;
%countpsi = 0;
countphih = 0;
countphig = 0;

%% initialize for storeage
stheta = zeros(nloop - burnin,6);  %[muh mug phih phig sigh2 sigg2]
%spsi = zeros(nloop - burnin,1); 
stau = zeros(nloop - burnin,T); 
sh = zeros(nloop - burnin,T);
sg = zeros(nloop - burnin,T);

%% compute a few things outside the loop
newnuh = T/2 + nuh0;
newnug = T/2 + nug0;
%psipri = @(x) -log(normpdf(x,0,sqrt(1/invVpsi))/(normcdf(sqrt(invVpsi))-normcdf(-sqrt(invVpsi))));

rand('state', sum(100*clock) ); randn('state', sum(200*clock) );

% pre-allocate
h_fore       = NaN(nloop-burnin,horizmax);
g_fore       = NaN(nloop-burnin,horizmax);
u_fore       = NaN(nloop-burnin,horizmax);
eps_y_fore   = NaN(nloop-burnin,horizmax);
y_fore       = NaN(nloop-burnin,horizmax);
tau_fore     = NaN(nloop-burnin,horizmax);
eps_tau_fore = NaN(nloop-burnin,horizmax);

mean_y_fore  = NaN(nloop-burnin,horizmax);
vari_y_fore  = NaN(nloop-burnin,horizmax);

for loop = 1:nloop

    %% sample tau    
    invS_tau = sparse(1:T,1:T,exp(-g));
    invOmega_tau = H'*invS_tau*H;  
    invS_y = sparse(1:T,1:T,exp(-h));
    invD_tautilde = invS_y + Hpsi'*invOmega_tau*Hpsi;
    tauhat = invD_tautilde\(invS_y*(Hpsi\y));
    tautilde = tauhat + chol(invD_tautilde)\randn(T,1);
    tau = Hpsi*tautilde;
    
    %% sample h
    Ystar = log((Hpsi\(y-tau)).^2 + .0001);
    h = SV(Ystar,h,phih,sigh2,muh);
   
    %% sample g
    Ystar = log([tau(1); tau(2:end)-tau(1:end-1)].^2 + .0001);
    g = SV(Ystar,g,phig,sigg2,mug);   
   
    %% sample sigh2
    errh = [(h(1)-muh)*sqrt(1-phih^2);  h(2:end)-phih*h(1:end-1)-muh*(1-phih)];
    newSh = Sh0 + sum(errh.^2)/2;
    sigh2 = 1/gamrnd(newnuh, 1./newSh);     
    
    %% sample sigg2
    errg = [(g(1)-mug)*sqrt(1-phig^2);  g(2:end)-phig*g(1:end-1)-mug*(1-phig)];
    newSg = Sg0 + sum(errg.^2)/2;
    sigg2 = 1/gamrnd(newnug, 1./newSg);

    %% sample psi    
%     fpsi = @(x) fMA1(x,y-tau,h) + psipri(x);
%     [psi flag psihat invDpsic] = sample_psi(psi,fpsi,loop,invDpsic,options);
%     Hpsi = speye(T) + sparse(2:T,1:(T-1),psi*ones(1,T-1),T,T); 
%     countpsi = countpsi + flag;
    
    %% sample phih
%     Xphih = h(1:end-1)-muh;
%     yphih = h(2:end) - muh;
%     Dphih = 1/(invVphih + Xphih'*Xphih/sigh2);
%     phihhat = Dphih*(invVphih*phih0 + Xphih'*yphih/sigh2);
%     phihc = phihhat + sqrt(Dphih)*randn;
%     fh = @(x) -.5*log(sigh2./(1-x.^2))-.5*(1-x.^2)/sigh2*(h(1)-muh)^2;
%     if abs(phihc)<.9999
%         alp = exp(fh(phihc)-fh(phih));
%         if alp>rand
%             phih = phihc;
%             countphih = countphih+1;
%         end
%     end 
    phih = 1;
    
    %% sample phig
%     Xphig = g(1:end-1)-mug;
%     yphig = g(2:end) - mug;
%     Dphig = 1/(invVphig + Xphig'*Xphig/sigg2);
%     phighat = Dphig*(invVphig*phig0 + Xphig'*yphig/sigg2);
%     phigc = phighat + sqrt(Dphig)*randn;
%     fg = @(x) -.5*log(sigg2./(1-x.^2))-.5*(1-x.^2)/sigg2*(g(1)-mug)^2;
%     if abs(phigc)<.9999
%         alp = exp(fg(phigc)-fg(phig));
%         if alp>rand
%             phig = phigc;
%             countphig = countphig+1;
%         end
%     end 
    phig = 1;
    
    %% sample muh    
%     Dmuh = 1/(invVmuh + ((T-1)*(1-phih)^2 + (1-phih^2))/sigh2);
%     muhhat = Dmuh*(invVmuh*muh0 + (1-phih^2)/sigh2*h(1) + (1-phih)/sigh2*sum(h(2:end)-phih*h(1:end-1)));
%     muh = muhhat + sqrt(Dmuh)*randn;   
    muh = 0;
    
    %% sample mug
%     Dmug = 1/(invVmug + ((T-1)*(1-phig)^2 + (1-phig^2))/sigg2);
%     mughat = Dmug*(invVmug*mug0 + (1-phig^2)/sigg2*g(1) + (1-phig)/sigg2*sum(g(2:end)-phig*g(1:end-1)));
%     mug = mughat + sqrt(Dmug)*randn;   
    mug = 0;
    
    if ( mod( loop, 5000 ) ==0 )
        disp(  [ num2str( loop ) ' loops... ' ] )
    end    
    
    if loop>burnin
        i = loop-burnin;
        stau(i,:) = tau';
        %spsi(i,:) = psi;
        sh(i,:) = h'; 
        sg(i,:) = g'; 
        stheta(i,:) = [muh mug phih phig sigh2 sigg2];     
        
        %forecast
        % horiz = 1
        resids = Hpsi\(y-tau);       
        h_fore(i,1) = muh + phih * ( h(end) - muh ) + sqrt(sigh2) * randn(1,1);
        g_fore(i,1) = mug + phig * ( g(end) - mug ) + sqrt(sigg2) * randn(1,1);
        
        u_fore(i,1)       = sqrt(exp(h_fore(i,1))) * randn(1,1);
        eps_tau_fore(i,1) = sqrt(exp(g_fore(i,1))) * randn(1,1);
        
        %eps_y_fore(i,1) = u_fore(i,1) + psi * resids(end); 
        eps_y_fore(i,1) = u_fore(i,1) +  0  * resids(end); 
        tau_fore(i,1) = tau(end) + eps_tau_fore(i,1);      
        
        y_fore(i,1) = tau_fore(i,1) + eps_y_fore(i,1);
        
        % calculate mean and variance forecast of each draw. These will be
        % used to determine the log score, since the forecast density will
        % be constructed as the linear pool of the normal densities of each
        % draw, where the normal density of draw i has the mean mean_y_fore(i,1) and the
        % variance vari_y_fore(i,1)
        mean_y_fore(i,1) = tau(end); % could have used                 = tau_fore(i,1)    and
        vari_y_fore(i,1)  =  exp(h_fore(i,1)) + exp(g_fore(i,1))  ;%   = exp(h_fore(i,1)) instead
        
        % horiz = 2,...
        for horiz = 2:horizmax
            h_fore(i,horiz) = muh + phih * ( h_fore(i,horiz-1) - muh ) + sqrt(sigh2) * randn(1,1);
            g_fore(i,horiz) = mug + phig * ( g_fore(i,horiz-1) - mug ) + sqrt(sigg2) * randn(1,1);
            
            u_fore(i,horiz)       = sqrt(exp(h_fore(i,horiz))) * randn(1,1);
            eps_tau_fore(i,horiz) = sqrt(exp(g_fore(i,horiz))) * randn(1,1);
            
            %eps_y_fore(i,horiz) = u_fore(i,horiz) + psi * u_fore(i,horiz-1);
            eps_y_fore(i,horiz) = u_fore(i,horiz)    + 0 * u_fore(i,horiz-1);
            tau_fore(i,horiz) = tau_fore(i,horiz-1) + eps_tau_fore(i,horiz);
            
            y_fore(i,horiz) = tau_fore(i,horiz) + eps_y_fore(i,horiz);
            
            mean_y_fore(i,horiz)  =  tau(end); 
            vari_y_fore(i,horiz)  =  exp(h_fore(i,horiz)) + sum(exp(g_fore(i,1:horiz)));
            
        end
        
    end    
end

disp( ['MCMC takes '  num2str( etime( clock, start_time) ) ' seconds' ] );
disp(' ' );

% tauhat = mean(stau)';
% tauCI = quantile(stau,[.05 .95])';
% hhat = mean(exp(sh/2))'; 
% hCI = quantile(exp(sh/2),[.05 .95])';
% ghat = mean(exp(sg/2))'; 
% gCI = quantile(exp(sg/2),[.05 .95])';
% 
% figure(1);
% hold on 
%     plotCI(tid,tauCI(:,1),tauCI(:,2));
%     plot(tid,tauhat); 
%     title('\tau_t'); xlim([tid(1)-.5 tid(end)+.5]); box off;
% hold off 
% 
% figure(2);
% subplot(2,1,1); 
% hold on 
%     plotCI(tid,hCI(:,1),hCI(:,2));
%     plot(tid,hhat); 
%     title('exp(h_t/2)'); xlim([tid(1)-.5 tid(end)+.5]); box off;
% hold off
% subplot(2,1,2); 
% hold on 
%     plotCI(tid,gCI(:,1),gCI(:,2));
%     plot(tid,ghat); 
%     title('exp(g_t/2)'); xlim([tid(1)-.5 tid(end)+.5]); box off;
% hold off

%figure; 
%hist(spsi,50);
%title('\psi_1'); box off;



