function [tau_draws,tau_f_draws,sigma_dtau_draws,sigma_eps_draws,g_eps_draws,g_dtau_draws,scl_eps_draws,sigmatotal_eps_draws] = ucsv_outlier(y,n_burnin,n_draws,g_eps_prior,g_dtau_prior,scl_eps_vec,prob_scl_eps_vec)
 % Model and notation
  %  y(t) = tau(t) + eps(t)
  %    tau(t) = tau(t-1) + dtau(t)
  %    eps and d_tau(t) are white noises with
  %     time varying standard deviations sigma_eps(t) and sigma_dtau(t)
  %     The log-volatilities follow random walk process with innovation standard deviations given by the parameter "g"
  %     The prior for g is given in the input matrices g_xxx_prior.  The first column are a grid of values for g and the 
  %     the second column are the probabilities
  %     This version of the program allows for fat tails in the distribution of eps by using a scale mixture of normals
  %      The scales (standard deviations) are given in scl_eps_vec and their prior probabilities are given in prob_scl_eps_vec
  %
  %     Note: to impose a particular value of g, use a point-mass prior  
  %     The output below are the draws of these variables
  %      tau_draws are the draws of tau
  %      tau_f_draws are "filtered" draws of tau -- but using the full-
  %        sample draws of sigma_eps and sigma_dtau
  %      sigma_eps_draws are the draws of sigma_eps
  %      sigma_dtau_draws are the draws of sigma_dtau
  %      scl_eps_draw are the draws of the scale of epsilon (normalized by volatility)
  % Mixture parameters for log chisquared (1) approximation 
  % r_p are probabilities
  % r_m are means
  % r_s are standard deviations
  
  %{
  % 2-component mixture used in Stock-Watson (2007)
  r_p = [0.086 (1-0.086)]';
  r_m = [-7.472 -0.698]';
  r_s = [1.411 1.411]';
  %}
  
  %{
  % 7-component mixture from Kim, Shephard, and Chib ReStud (1998);
  r_p = [0.00730 0.10556 0.00002 0.04395 0.34001 0.24566 0.25750]';
  r_m = [-10.12999 -3.97281 -8.56686 2.77786 0.61942 1.79518 -1.08819]';
  r_v = [5.79596 2.61369 5.17950 0.16735 0.64009 0.34023 1.26261]';
  r_m = r_m - 1.2704;
  r_s = sqrt(r_v);
  %}
  
  % 10-component mixture from Omori, Chib, Shephard, and Nakajima JOE (2007)
  r_p = [0.00609 0.04775 0.13057 0.20674 0.22715 0.18842 0.12047 0.05591 0.01575 0.00115]';
  r_m = [1.92677 1.34744 0.73504 0.02266 -0.85173 -1.97278 -3.46788 -5.55246 -8.68384 -14.65000]';
  r_v = [0.11265 0.17788 0.26768 0.40611 0.62699 0.98583 1.57469 2.54498 4.16591 7.33342]'; 
  r_s = sqrt(r_v);
  
  T = size(y,1);   % Sample size
  
  % Scale y:  This scale matters because of the off-set parameters "c" in 
  % draw_lcs_indicators and draw_sigma;
  dy = y(2:end)-y(1:end-1);
  inan = isnan(dy);
  std_dy = std(dy(inan==0));
  scale_y = std_dy/5;
  yn = y/scale_y;
  
  % Matrices for saving results
  tau_draws = NaN*zeros(T,n_draws);
  tau_f_draws = NaN*zeros(T,n_draws);
  sigma_dtau_draws = NaN*zeros(T,n_draws);
  sigma_eps_draws = NaN*zeros(T,n_draws);
  g_eps_draws = NaN*zeros(n_draws,1);
  g_dtau_draws = NaN*zeros(n_draws,1);
  scl_eps_draws = NaN*zeros(T,n_draws);
  sigmatotal_eps_draws = NaN*zeros(T,n_draws);
  
  % -- Begin Calculations: sigma = 1 and scale_eps is unity;
  sigma_eps = ones(T,1);
  sigma_dtau = ones(T,1);
  scale_eps = ones(T,1);
  
  itmp = 0;
  tic;
  for idraw=1:n_burnin+n_draws;
  	% Step 1(a): draw tau, tau_f, dtau, and eps; 
  	sigma_eps_scl = sigma_eps.*scale_eps;	  % SD of eps, which is stochastic volatility times scale in mixture distribution 
    [tau_a,tau_f] = draw_tau(yn,sigma_dtau,sigma_eps_scl);
    dtau = tau_a(2:end)-tau_a(1:end-1);
    tau = tau_a(2:end);
    eps = yn - tau;
    eps_scaled = eps./scale_eps;  % Scaled version of epsilon
  
    % Step 1(b): Draw mixture indicators for log chi-squared(1)
    ind_eps = draw_lcs_indicators(eps_scaled,sigma_eps,r_p,r_m,r_s); 
    ind_dtau = draw_lcs_indicators(dtau,sigma_dtau,r_p,r_m,r_s);
    
    % Step 2(a): Draw G
    g_eps = draw_g(eps_scaled,g_eps_prior,ind_eps,r_m,r_s,1);
    g_dtau = draw_g(dtau,g_dtau_prior,ind_dtau,r_m,r_s,1);
    
    % Step 2(b): Draw Volatilities
    sigma_eps = draw_sigma(eps_scaled,g_eps,ind_eps,r_m,r_s,1);	
    sigma_dtau = draw_sigma(dtau,g_dtau,ind_dtau,r_m,r_s,1);	
    
    % Step 3: Draw Scale of epsilon
    
    scale_eps = draw_scale_eps(eps,sigma_eps,ind_eps,r_m,r_s,scl_eps_vec,prob_scl_eps_vec);
   
    % Save draws;
    if idraw > n_burnin;
      j = idraw- n_burnin;
      tau_draws(:,j) = tau;
      tau_f_draws(:,j) = tau_f;
      sigma_dtau_draws(:,j) = sigma_dtau;
      sigma_eps_draws(:,j) = sigma_eps;
      g_eps_draws(j) = g_eps;
      g_dtau_draws(j) = g_dtau;
      scl_eps_draws(:,j) = scale_eps;
      sigmatotal_eps_draws(:,j) = sigma_eps.*scale_eps;
    end;
    
    itmp=itmp+1;
    if itmp == 1000;
      itmp=0;
      idraw;
      toc
      tic;
    end;
    
  end;
  
  % Rescale so that things are in units of y
  tau_draws = tau_draws*scale_y;
  tau_f_draws = tau_f_draws*scale_y;
  sigma_dtau_draws = sigma_dtau_draws*scale_y;
  sigma_eps_draws = sigma_eps_draws*scale_y;
  sigmatotal_eps_draws = sigmatotal_eps_draws*scale_y;

end

