function [st_f, theta_f]=volatility_LP(vardata,p,c,T0,LP)

[Y,X,Y_initial]=SUR(vardata,p,c);

[T,n]=size(Y);

TT=T0-p; % last pre-COVID obs - p

    if LP==0
        st_f = ones(size(Y,1),1);
        theta_f=[];

    elseif LP==1
    nst=3;
    
    theta0 = ones(nst,1);	% initialize parameters vectors
    theta0(1,1) =15;
    theta0(2,1) =15;
    theta0(3,1) =0.7;
    
    f = @(theta) LP_Covid(theta,Y,TT,n,p, X);	% define log-likelihood
    
    theta_f = estTheta(f,theta0,nst);  
    
    st_f = build_st(theta_f,T,TT); % vector S
    
    elseif LP==2
        nst = T-TT; % nst = n. of parameters to estimate (one parameter for each post-covid obs)
                    
        theta0 = ones(nst,1);	% initialize parameters vector
        f = @(theta) LP_Covid(theta,Y,TT,n,p, X);	% define log-likelihood

    
    theta_f = estTheta(f,theta0,nst); % Estimate Covid volatility parameters
    
    st_f = build_st(theta_f,T,TT); % Create S vector
    end

end
