function [THETA,THETA_sum,elapsed_time,keeps_b,Ypred,Ypred_sum] = estim_CEVe(Y,m,H,hor,THETA,flag_beta)
% ESTIM_CEVe estimates the parameters of the CEV model as described in
%   Eraker (2001), using MCMC methods. Potentially, missing data is allowed.

% REQUIRED INPUT
% 'Y' is the (n X 1) vector with data.
% 'Dt = 1/m' is increment in data, thus m-1 is number of missing values.
%   NOTE: this program only deals with m = 1.
% 'H' can be an scalar with the total number of iterations,
%   or a 1x2 vector [H0,Htot], where H0 are the 'burn-in' iterations, and
%   the Htot the total number of iterations in the gibbs sampler and MH algorithm.
%   Thus, the posterior means of the parameters are calculated with 
%   Hr := Htot-H0
%   observations. If H is a scalar H0 = ceiling(H/4).
% OPTIONAL INPUT
% 'hor' is the forcasting horizont. The default is hor = 0.
% 'THETA' is a (1 x cT) vector of initial values for parameters, where 
%   cT = 3,4, depending on the model. The default is a (1x4) vector.
% 'flag_beta' is a logical 0/1 variable that enables (1) or disablabes (0) the
%   estimation of beta. It is also used to define the model estimated.
%   The default depends on size of THETA. If THETA is (1x3) flag_beta = 0; 
%   o.w. the default is flag_beta = 1. 
%
% OUTPUT VARIABLES
% 'THETA' is the (Hr x cT) matrix of draws from parameter posterior distributions.
% 'THETA_sum' is a (3 x cT) summary matrix: 1st row means, 2nd row standard errors,
%   and 3rd logical vector specifying whether a paramter was estimated (1) or not (0).
% 'elapsed_time' records the time the computation time of the MCMC. 
% 'keeps_b' is a (1 x Htot) vector with the counter of draws from beta kept.
%   It gives info about proficiency of MH.
% 'Ypred' is a (hor x Hc) matrix with predictive posterior draws. 
% 'Ypred_sum' is a (2 x hor) matrix with the corresponding
%   predictive posterior means and standard deviations from the above matrix.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CEV model is                                                            %
% r(t)-r(t-1) = mu(r(t-1),theta)*Dt + sigma(r(t-1),theta)*sqrt(Dt)*e^y(t) %
% where:                                                                  %
% mu(r_t,theta) = thetaR+kappaR*r_t,                                      %
% sigma(r_t,theta) = sigmaR*r_t^beta,                                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin < 4
    hor = 0;
end
if nargin < 5
    THETA = [0.01,-0.02,0.5,0.5];
end
if nargin < 6
    if length(THETA) == 3
        flag_beta = 0;
    else
        flag_beta = 1;
    end
end

temp0 = Y;
Ti=size(temp0,1);
m = 1;
Dt=1/m;
ind1=(1:m:m*Ti);
n=ind1(Ti);
ind0=1:n;

if m==1
    Y=temp0';
else
    Y=zeros(size(temp0,2),ind1(Ti));
    indt=(ind1+m-1)/m;
    Y(:,ind1)=temp0(indt,:)';
    ind2 = ind_f(ind1,m,Ti);
    for j = 1:(m-1)
        Y(:,j+ind1(1:(Ti-1))) = (Y(:,j+ind1(1:Ti-1)-1)+Y(:,ind1(1:Ti-1)+m))/2;
    end
end

t1 = datevec(now);  % time of start of computation
[THETA,THETA_sum,keeps_b,Ypred,Ypred_sum] = estim_CEVf(Y,m,H,hor,THETA,flag_beta);
t2 = datevec(now);  % time of end of computation
[elapsed_time] = took(t1,t2);


%%%%%%%%%%%%%%%%%%%%
% MAIN SUBFUNCTION %
%%%%%%%%%%%%%%%%%%%%

function [THETA,THETA_sum,keeps_b,Ypred,Ypred_sum] = estim_CEVf(Y,m,H,hor,THETA,flag_beta)
% This function gives all models, depending on several variables.

if length(H)==1
    H0 = ceil(H/4);
else
    H0 = H(1);
    H = H(2);
end

Ti=size(Y,2);
m = 1;
Dt=1/m;
ind1=(1:m:m*Ti);
n=ind1(Ti);
ind0=1:n;

% initialization of parameters   
thetaR = THETA(1);
kappaR = THETA(2);
sigmaR = THETA(3);
if length(THETA) == 3
    beta = 0;
else
    beta = THETA(4);
end

keeps_b =[];
Ypred = [];
Ypred_h = zeros(1,hor+1);  % if Y is multirow, change to (d1,hor+1)
Ypred_h(:,1) = Y(:,n);

if m==1  
    % No augmented data, therefore skip to step 2 in pp 182    
    for h = 2:H
        [tmp1,tmp2,tmp3] = paramY(n,Y,beta(h-1),sigmaR(h-1));
        thetaR = [thetaR; tmp1];
        kappaR = [kappaR; tmp2];
        sigmaR = [sigmaR; tmp3];

        % Metropolis-Hastings for beta
        if flag_beta == 1
            newVal = q_beta_r(beta(h-1),Y,ind1(Ti),sigmaR(h),thetaR(h),kappaR(h),Dt);
            temp5 = lp_beta(newVal,Y,ind1(Ti),sigmaR(h),thetaR(h),kappaR(h),Dt)-lp_beta(beta(h-1),Y,ind1(Ti),sigmaR(h),thetaR(h),kappaR(h),Dt);
            alpha_b = exp(temp5)*q_beta(beta(h-1),beta(h-1),Y,ind1(Ti),sigmaR(h),thetaR(h),kappaR(h),Dt)/q_beta(newVal,beta(h-1),Y,ind1(Ti),sigmaR(h),thetaR(h),kappaR(h),Dt);
            alpha_b = min1(alpha_b);
            keep = binornd(1,alpha_b);  % indicates if you keep value
            beta = cat(1,beta,beta(h-1)+(newVal-beta(h-1))*keep);
            keeps_b = cat(2,keeps_b,keep);  % counter of values kept
        else
            beta = [beta; beta(h-1)];
        end
        
        % Prediction
        if hor > 0
            for k = 2:(hor+1)
                m1 = Ypred_h(:,k-1)+(thetaR(h)+kappaR(h)*Ypred_h(:,k-1))*Dt; 
                sd1 = sigmaR(h)*(abs(Ypred_h(:,k-1))^beta(h))*sqrt(Dt);
                temp6 = normrnd(m1,sd1);
                if flag_beta == 1
                    while temp6 <= 0
                        temp6 = normrnd(m1,sd1);
                    end
                end
                Ypred_h(:,k) = temp6;
            end
            Ypred = cat(1,Ypred,Ypred_h(2:hor+1));
        end

    end
else
    error('Missing values not allowed: "m" must be equal to 1')
end

if hor > 0
    Ypred_sum = cat(1,mean(Ypred(H0:(H-1),:)),std(Ypred(H0:(H-1),:)));
else
    Ypred_sum = [];
end


%THETA = cat(2,thetaR,kappaR,sigmaR,beta);
THETA = cat(2,1000*thetaR,kappaR,sigmaR,beta);
THETA = THETA((H0+1):H,:);
estim = [1 1 1 flag_beta];
THETA_sum = cat(1,mean(THETA),std(THETA),estim);
%


%%%%%%%%%%%%%%%%
% subfunctions %
%%%%%%%%%%%%%%%%

function Dy = diff_s(y)
% DIFF_S first difference of time series vector
% DIFF_S(y), where y is a mxn matrix and n is no. of observations
n = size(y,1);
temp1 = y(1:(n-1),:);
temp2 = y(2:n,:);
Dy = temp2 - temp1;
%

function [beta0,beta1,sig] = paramY(n,Y,beta,sigma)
% PARAMy samples from posterior distribution of mean equation of returns
% The result is simple because thetaK,kappa ~ N(0,I)
if beta == 0
    temp1 = (diff_s(Y')');
    temp2 = cat(1,ones(1,(n-1)),Y(:,1:(n-1)));
else
    temp1 = (diff_s(Y')')./(abs(Y(:,1:(n-1))).^beta);
    temp2 = cat(1,Y(:,1:(n-1)).^(-beta),Y(:,1:(n-1)).^(1-beta));
end
y = temp1';
X = temp2';
Sig = inv(X'*X);
theta_bar = Sig*X'*y;
res = y-X*theta_bar; 
s_bar2 = (1/n)*res'*res;

temp3 = mvnrnd(theta_bar,sigma^2*Sig,1); 
beta0 = temp3(1);
beta1 = temp3(2);
sig = sqrt(1/gamrnd(n/2-1,1/(n*s_bar2/2)));
%

function lpbet = lp_beta(beta,Y,n,sigmaR,thetaR,kappaR,Dt)
% LP_BETA log-posterior of beta given data: -beta*ln(x)-1/2*c*x^(-2*beta)
temp1 = diff_s(Y(1,:)')';
temp2 = Y(1,1:(n-1));
cons = ((temp1-(thetaR+kappaR.*temp2).*Dt).^2)./((sigmaR.^2)*Dt); % i removed 0.5*
[m1,n1] = size(temp1);
if m1 < n1
    lpbet = sum(-beta.*log(temp2)-(1/2)*cons.*temp2.^(-2*beta),2);
else
    lpbet = sum(-beta.*log(temp2)-(1/2)*cons.*temp2.^(-2*beta),1);
end
%

function [lbet,lbet2] = lp_dbeta(beta,Y,n,sigmaR,thetaR,kappaR,Dt)
% LP_BETA first derivative of the log-posterior: -ln(x)+c*x^(-2*beta)*ln(x)
% LP_BETA2 second derivative of the log-posterior: c*x^(-2*beta)*(ln(x))^2
temp1 = diff_s(Y')';
temp2 = Y(1:(n-1));
cons = ((temp1-(thetaR+kappaR.*temp2).*Dt).^2)./((sigmaR.^2)*Dt); % changed from 0.25*Z on 08/20/2003
[m1,n1] = size(temp1);
if m1 < n1
    lbet = sum(-log(temp2)+cons.*(temp2.^(-2*beta)).*log(temp2),2);
    lbet2 = -2*sum(cons.*temp2.^(-2*beta).*(log(temp2)).^2,2);
else
    lbet = sum(-log(temp2)+cons.*(temp2.^(-2*beta)).*log(temp2),1);
    lbet2 = -2*sum(cons.*temp2.^(-2*beta).*(log(temp2)).^2,1);
end
%

function val = q_beta(x,beta,Y,n,sigma,thetaR,kappaR,Dt)
% Q_BETA evaluates x at beta posterior proposal density
m = mean_b(beta,Y,n,sigma,thetaR,kappaR,Dt);
s2 = sigma2_b(beta,Y,n,sigma,thetaR,kappaR,Dt);
val = normpdf(x,m,sqrt(s2));
%

function ran_b = q_beta_r(beta,Y,n,sigma,thetaR,kappaR,Dt)
% Q_BETA_R generates a random number from beta posterior proposal
m = mean_b(beta,Y,n,sigma,thetaR,kappaR,Dt);
s2 = sigma2_b(beta,Y,n,sigma,thetaR,kappaR,Dt);
ran_b = normrnd(m,sqrt(s2));
%

function mbet = mean_b(beta,Y,n,sigma,thetaR,kappaR,Dt)
% MEAN_B mean of proposal for beta posterior
[a,b] = lp_dbeta(beta,Y,n,sigma,thetaR,kappaR,Dt);
mbet = beta-a/b;
%

function s2bet = sigma2_b(beta,Y,n,sigma,thetaR,kappaR,Dt)
% SIGMA2_B variance of proposal for beta posterior
[c,d] = lp_dbeta(beta,Y,n,sigma,thetaR,kappaR,Dt);
s2bet = -1/d;
%

function mx = min1(x)
%MIN1 returns the minimum between abs(x) and 1
if abs(x) < 1
    mx = abs(x);
else
    mx = 1;
end
%

function el_t = took(t1,t2)
% TOOK computes how many minutes or seconds a computation took
tdif = t2 - t1;
temp = tdif(3)*24*60*60 + tdif(4)*3600 + tdif(5)*60 +tdif(6);
if temp < 60
    th = 0;
    tm = 0;
    ts = temp;
elseif temp < 3600
    th = 0;
    tmp1 = temp/60;
    tm = floor(tmp1);
    ts = (60*(tmp1-tm));
else
    tmp1 = temp/3600;
    th = floor(tmp1);
    tmp2 = (60*(tmp1-th));
    tm = floor(tmp2);
    ts = (60*(tmp2-tm));
end
el_t = sprintf('%1.0f hr %1.0f min %1.2f sec',th,tm,ts);
%