function [tau]=ucsv(y);

%------------------------------------------------------------%
% This script provides an MCMC algorithm to estimate an      %
% univariate UC-SV model a la Stock and Watson (2007) using  %
% precision sampling methods as in Chan and Jeliazkov (2009) %
%                                                            %
%                      Written by Luis Uzeda (February-2016) %
%------------------------------------------------------------%

% Preamble
% clear all; clc;
% load 'PCE_US.csv' %upload dataset 
% dta=PCE_US; %[PCE-Inflation, PCE-Inflation-Goods, PCE-Inflation-Services, Goods Share, Services Share]
% ts=1;       %select inflation series; 1=Aggregate Inflation, 2=Goods Inflation, 3=Services Inflation 
% y=dta(:,ts);

yOrig=y;
nanIndex=isnan(y);
y=y(~nanIndex);

nloop=8000;
burnin=2000;
T=length(y);

%allocating storage
store_eta=zeros(T,(nloop-burnin));  
store_tau=zeros(T,(nloop-burnin));
store_heps=zeros(T,(nloop-burnin));
store_heta=zeros(T,(nloop-burnin));
store_theta=zeros(nloop-burnin,2); %[sig2eps sig2eta]

%% =========================Priors===========================%

%----sig2varxieps, sig2varxieta-----% 
%if using inverse-Gamma prior for variance parameters
nusig2varxiepspr=15; Ssig2varxiepspr=(.17)^2*(nusig2varxiepspr-1); 
nusig2varxietapr=15; Ssig2varxietapr=(.17)^2*(nusig2varxietapr-1);

%------z0 (initial conditions (tau))-----%
tau0=.99*y(1);
z0=tau0;

%-----set variance to initialize log-volatilities-----% 
sig2varxieps0=10^7; 
sig2varxieta0=10^7;

%===============Initial values for the Markov Chain==============%
invSigeps=speye(T)*.1;
invSigeta=speye(T)*.1;
sig2varxieps=.01;
sig2varxieta=.01;
sig2varxiepsold=.224^2;
sig2varxietaold=.224^2;
heps=.5*ones(T,1);
heta=.5*ones(T,1);

%=======Compute a few other things outside the MCMC loop========%
H=speye(T)+sparse(2:T,1:T-1,-1*ones(T-1,1),T,T);
iota0=[1;zeros(T-1,1)];

%% begin MCMC
start_time = clock;
rand('state', sum(100*clock) ); randn('state', sum(200*clock) );
disp(' ' );
for iter=1:nloop
%% (1)-Disturbance Smoothing

%sample eta (vector of innovations to trend inflation)
A=H;
X0=A\iota0;
X1=speye(T);
ytemp=H*y;
ytil=A\ytemp;
yetatil=ytil-X0*z0;
invDetapr=A'*invSigeta*A;
invDetatil=X1'*invSigeps*X1+invDetapr;
C=chol(invDetatil,'lower');
etatilnew=C'\(C\(X1'*invSigeps*yetatil));
etatil=etatilnew+C'\randn(T,1);
eta=A*etatil;

%% Sampling log-volatilities 

%---------espilon (innovations to transitory component)--------%
%sample s from a 7-point discrete distriburion (Kim et al (1998))
eps=ytil-X0*z0-X1*etatil;  
ytilh = log(eps.^2 + .0001 );
% normal mixture 
pmix = [0.0073 .10556 .00002 .04395 .34001 .24566 .2575];
mi = [-10.12999 -3.97281 -8.56686 2.77786 .61942 1.79518 -1.08819] - 1.2704;  
sigi = [5.79596 2.61369 5.17950 .16735 .64009 .34023 1.26261];
sqrtsigi = sqrt(sigi);
temprand = rand(T,1);
liktemp=normpdf(repmat(ytilh,1,7),repmat(heps,1,7)+repmat(mi,T,1),repmat(sqrtsigi,T,1));
q = repmat(pmix,T,1).*liktemp; %prob of drawing s_t|y
q = q./repmat(sum(q,2),1,7); %normalizing prob of drawing s|y
s1 = 7 - sum(repmat(temprand,1,7)<cumsum(q,2),2)+1;

% sample h_eps
invSigvarxieps = spdiags([1/(sig2varxieps0); 1/(sig2varxieps*1)*ones(T-1,1)],0,T,T);
alph_mix = mi(s1)';
ytilhstar=ytilh-alph_mix;
invSigs = spdiags(1./sigi(s1)',0,T,T);
invSighnew = invSigs + H'*invSigvarxieps*H;
C = chol(invSighnew,'lower');
hepsnew = C'\(C\(invSigs*(ytilhstar)));
heps = hepsnew + C'\randn(T,1);

%covariance matrix-epsilon
invSigeps=spdiags(exp(-heps),0,T,T);

%------------eta (innovations to trend)------------%
%sample s from a 7-point discrete distriburion (Kim et al (1998))
etastr = log(eta.^2 + .0001 );
% normal mixture 
pmix = [0.0073 .10556 .00002 .04395 .34001 .24566 .2575];
mi = [-10.12999 -3.97281 -8.56686 2.77786 .61942 1.79518 -1.08819] - 1.2704;  
sigi = [5.79596 2.61369 5.17950 .16735 .64009 .34023 1.26261];
sqrtsigi = sqrt(sigi);
temprand = rand(T,1);
q = repmat(pmix,T,1).*normpdf(repmat(etastr,1,7),repmat(heta,1,7)+repmat(mi,T,1),repmat(sqrtsigi,T,1)); %prob of drawing s_t|y
q = q./repmat(sum(q,2),1,7); %normalizing prob of drawing s|y
s2 = 7 - sum(repmat(temprand,1,7)<cumsum(q,2),2)+1;  

% sample h_eta
invSigvarxieta = spdiags([1/sig2varxieta0; 1/sig2varxieta*ones(T-1,1)],0,T,T);
alph_mix = mi(s2)';
eta2str=etastr-alph_mix;
invSigs = spdiags(1./sigi(s2)',0,T,T);
invSighetanew = invSigs + H'*invSigvarxieta*H;
C = chol(invSighetanew,'lower');
hetanew = C'\(C\(invSigs*(eta2str)));
heta = hetanew + C'\randn(T,1);

%covariance matrix-eta
invSigeta=spdiags(exp(-heta),0,T,T);


%% Parameter Sampling

%sample sig2varxieps
nusig2varxieps=(T-1)/2+nusig2varxiepspr;
Ssig2varxieps=sum((heps(2:end)-heps(1:end-1)).^2)/2+Ssig2varxiepspr;   
sig2varxieps=1/gamrnd(nusig2varxieps,1/Ssig2varxieps);

%sample sig2varxieta
nusig2varxieta=(T-1)/2+nusig2varxietapr;
Ssig2varxieta=sum((heta(2:end)-heta(1:end-1)).^2)/2+Ssig2varxietapr;  
sig2varxieta=1/gamrnd(nusig2varxieta,1/Ssig2varxieta);    


%% Recovering Trend Inflation (tau=y-epsilon)
tau=y-eps;

%% Begin storage  
      if iter>burnin 
          i = iter-burnin;            
          store_eta(:,i)=eta; 
          store_tau(:,i)=tau;
          store_heps(:,i)=heps;
          store_heta(:,i)=heta;
          store_theta(i,:) = [sig2varxieps sig2varxieta];
      end
    if ( mod( iter, 2000 ) ==0 )
        disp(  [ num2str( iter ) ' iterations... ' ] )
    end
end    
disp( ['MCMC takes '  num2str( etime( clock, start_time) ) ' seconds' ] );
disp(' ' );   

% 
% %% Printing Results
% 
% %------Plotting Measures of Trend-Inflation------%
% figure;
% t=1947.25:.25:2017;
% g1=20;g2=20;g3=20;
% tau=median(store_tau'); %#ok
% q05=quantile(store_tau',.05);
% q95=quantile(store_tau',.95);
% hold 'on';
% %Plot shaded areas for NBER recession periods
% area([1948.75,1949.75],[g1,g2],-g3,'FaceColor',.8*[1,1,1],'EdgeColor', 'none');
% area([1953.25,1954.25],[g1,g2],-g3,'FaceColor',.8*[1,1,1],'EdgeColor', 'none');
% area([1957.5,1958.25],[g1,g2],-g3,'FaceColor',.8*[1,1,1],'EdgeColor', 'none');
% area([1960.25,1961],[g1,g2],-g3,'FaceColor',.8*[1,1,1],'EdgeColor', 'none');
% area([1969.75,1970.75],[g1,g2],-g3,'FaceColor',.8*[1,1,1],'EdgeColor', 'none');
% area([1973.75,1975],[g1,g2],-g3,'FaceColor',.8*[1,1,1],'EdgeColor', 'none');
% area([1980,1980.5],[g1,g2],-g3,'FaceColor',.8*[1,1,1],'EdgeColor', 'none');
% area([1981.5,1982.75],[g1,g2],-g3,'FaceColor',.8*[1,1,1],'EdgeColor', 'none');
% area([1990.5,1991],[g1,g2],-g3,'FaceColor',.8*[1,1,1],'EdgeColor', 'none');
% area([2000,2000.75],[g1,g2],-g3,'FaceColor',.8*[1,1,1],'EdgeColor', 'none');
% area([2007.75,2009.25],[g1,g2],-g3,'FaceColor',.8*[1,1,1],'EdgeColor', 'none');
% infl_agg=plot(t,dta(:,1),'g','LineWidth',2);
% trendinfl_agg=plot(t,tau,'b','LineWidth',2);
% q5=plot(t,q05,'--r','LineWidth',1);
% q95=plot(t,q95,'--r','LineWidth',1);
% plot(get(gca,'Xlim'),[0 0],'k-');
% lgnd=legend([infl_agg,trendinfl_agg,q5,q95],'Headline Inflation', 'Trend Inflation','5%-tile','95%-tile','Location','NorthEast');
% set(lgnd, 'Box', 'off');
% hold 'off';
% 
% %------Plotting Measures of Volatility------%
% figure;
% t=1947.25:.25:2017;
% g1=6;g2=6;g3=-0;
% vol_trend=exp(median(store_heta')./2); %#ok
% vol_trans=exp(median(store_heps')./2); %#ok
% hold 'on';
% %Plot shaded areas for NBER recession periods
% area([1948.75,1949.75],[g1,g2],-g3,'FaceColor',.8*[1,1,1],'EdgeColor', 'none');
% area([1953.25,1954.25],[g1,g2],-g3,'FaceColor',.8*[1,1,1],'EdgeColor', 'none');
% area([1957.5,1958.25],[g1,g2],-g3,'FaceColor',.8*[1,1,1],'EdgeColor', 'none');
% area([1960.25,1961],[g1,g2],-g3,'FaceColor',.8*[1,1,1],'EdgeColor', 'none');
% area([1969.75,1970.75],[g1,g2],-g3,'FaceColor',.8*[1,1,1],'EdgeColor', 'none');
% area([1973.75,1975],[g1,g2],-g3,'FaceColor',.8*[1,1,1],'EdgeColor', 'none');
% area([1980,1980.5],[g1,g2],-g3,'FaceColor',.8*[1,1,1],'EdgeColor', 'none');
% area([1981.5,1982.75],[g1,g2],-g3,'FaceColor',.8*[1,1,1],'EdgeColor', 'none');
% area([1990.5,1991],[g1,g2],-g3,'FaceColor',.8*[1,1,1],'EdgeColor', 'none');
% area([2000,2000.75],[g1,g2],-g3,'FaceColor',.8*[1,1,1],'EdgeColor', 'none');
% area([2007.75,2009.25],[g1,g2],-g3,'FaceColor',.8*[1,1,1],'EdgeColor', 'none');
% vol_trend=plot(t,vol_trend,'b','LineWidth',2);
% vol_trans=plot(t,vol_trans,'--r','LineWidth',2);
% set(gca,'Xlim',[1947.25,2016.75]);
% legend([vol_trend,vol_trans], 'Trend Inflation', 'Transitory Inflation', 'Location','NorthWest')
% legend('boxoff');
% title('Time Varying Volatility')
% hold 'off';
% 
% 
% 
% 
% 



