mrs_est.m

Modeling and Forecasting Electricity Loads and Prices: A Statistical Approach" by Rafa&#322 Weron, p

                    logL,   nan,    nan,       nan,    nan,     nan,    nan];
end

%=========================================================================
% mrs_Moments
%=========================================================================

function [sMean, sVar] = mrs_Moments(Model, sParam)

if strcmp(Model, 'G') == 1
    %Gaussian distribution:
    sC = sParam(2);
    sSigma2 = sParam(3);
    sMean = sC;
    sVar = sSigma2;
elseif strcmp(Model, 'LN') == 1
    %log-normal distribution
    sC = sParam(2);
    sSigma2 = sParam(3);
    sMean = exp(sC + sSigma2/2);
    sVar = exp(2*sC+sSigma2)*(exp(sSigma2)-1);
elseif strcmp(Model, 'P') == 1
    %Pareto distribution
    sA = sParam(4);
    sK = sParam(5);
    sMean = (sA*sK)/(sA-1);
    sVar = (sA*sK*sK)/((sA-1)*(sA-1)*(sA-2));
elseif strcmp(Model, 'W') == 1
    %Weibull distribution
    sA = sParam(4);
    sB = sParam(5);
    sMean = sA^(-1/sB)*gamma(1+1/sB);
    sVar = sA^(-2/sB)*(gamma(1+2/sB)-(gamma(1+1/sB))^2);
elseif strcmp(Model, 'MR') == 1
    sPhi = sParam(1);
    sC = sParam(2);
    sSigma2 = sParam(3);
    sMean =  sC / (1-sPhi);
    sVar = sSigma2 / (1-sPhi^2);
end

%=========================================================================
% mrs_EM_MLE
%=========================================================================

function newParam = mrs_EM_MLE(Data, Model, Ksi_tT, oldParam)
% MRS_EM_MLE ML estimation of regime processes parameters. Note that 
% Pareto distribution uses moment estimators.

sData = Data(2:end);
sKsi_tT = Ksi_tT(2:end, 2);

bData = Data(1:end);
bKsi_tT = Ksi_tT(2:end, 1);

% Spike regime 
switch Model
    case 'G'
        %Gaussian distribution
        sNewC = sum(sData.*sKsi_tT) / sum(sKsi_tT);
        sNewSigma2 = sum((sData-sNewC).^2 .* sKsi_tT) / sum(sKsi_tT);
        sNewParam = [nan, sNewC, sNewSigma2, nan, nan];
    case 'LN'
        %log-normal distribution
        sNewC = sum(log(sData).*sKsi_tT) / sum(sKsi_tT);
        sNewSigma2 = sum((log(sData)-sNewC).^2 .* sKsi_tT) / sum(sKsi_tT);
        sNewParam = [nan, sNewC, sNewSigma2, nan, nan];
    case 'P'
        %Pareto distribution
        sN = length(sData);
        sMean = sum(sData.*sKsi_tT) / sum(sKsi_tT);
        sMin = min(sData);
        % uncomment the next 2 lines to use ML estimators
        %  sNewA = (sMin-sN*sMean) / (sN*(sMin-sMean));
        %  sNewK = (sNewA*sN-1)*sMin/(sNewA*sN);
        % uncomment the next 2 lines to use moment estimators
        sNewK = sMin;
        sNewA = sum(sKsi_tT)/sum(sKsi_tT .* (log(sData)-log(sNewK)));
        sNewParam = [nan, nan, nan, sNewA, sNewK];       
    case 'W'
        %Weibull distribution
        sOldA = oldParam(2,4);
        sOldB = oldParam(2,5);
        sNewParam = fsolve(@weibullEM_MLE, [sOldA, sOldB], optimset('Display', 'off'), sData, sKsi_tT);
        sNewA = sNewParam(1);
        sNewB = sNewParam(2);
        sNewParam = [nan, nan, nan, sNewA, sNewB];
    case 'MR'
        %mean-reverting process
        sOldPhi = oldParam(2,1);
        sOldC = oldParam(2,2);
        sNewParam = fsolve(@mrEM_MLE, [sOldPhi, sOldC], optimset('Display', 'off'), bData, sKsi_tT);    
        sNewPhi = sNewParam(1);
        sNewC = sNewParam(2);
        sNewSigma2 = sum((sData - sNewC - sNewPhi*bData(1:end-1)).^2 .* sKsi_tT) / sum(sKsi_tT);
        sNewParam = [sNewPhi, sNewC, sNewSigma2, nan, nan];
end

% mean-reverting regime 
bOldPhi = oldParam(1,1);
bOldC = oldParam(1,2);
bNewParam = fsolve(@mrEM_MLE, [bOldPhi, bOldC], optimset('Display', 'off'), bData, bKsi_tT);
bNewPhi = bNewParam(1);
bNewC = bNewParam(2);
bNewSigma2 = sum((sData - bNewC - bNewPhi*bData(1:end-1)).^2 .* bKsi_tT) / sum(bKsi_tT);
bNewParam = [bNewPhi, bNewC, bNewSigma2, nan, nan];

newParam = [bNewParam; sNewParam];

%==========================================================================
% mrs_EM_MLE auxiliary functions 
%==========================================================================

function F = mrEM_MLE(mrParam, mrData, mrKsi_tT)
mrPhi = mrParam(1);
mrC = mrParam(2);
F = [sum((mrData(2:end)-mrC-mrPhi*mrData(1:end-1)).*mrKsi_tT);
    sum((mrData(2:end)-mrC-mrPhi*mrData(1:end-1)).*mrKsi_tT.*mrData(1:end-1))];

function F = weibullEM_MLE(wParam, wData, wKsi_tT)
wA = wParam(1);
wB = wParam(2);
F = [sum((1/wA - wData.^wB).*wKsi_tT);
    sum((1/wB + log(wData) - wA * log(wData) .* wData.^wB).*wKsi_tT)];

%==========================================================================
% mrs_EM_Smoother 
%==========================================================================

function [Ksi_tT, newP, newKsi_t1t_10, LogL] = mrs_EM_Smoother(Data, Model, oldParam, oldP, oldKsi_t1t_10)

N = length(Data);
    
% old base regime parameters
bOldPhi = oldParam(1, 1);
bOldC = oldParam(1, 2);
bOldSigma2 = oldParam(1, 3);

% Cond. dens. for the base regime (first column) and spike regime (second column)
switch Model
    case 'G'
        %Gaussian distribution
        sOldC = oldParam(2,2);
        sOldSigma2 = oldParam(2,3);
        Eta = [nan, nan; normpdf(Data(2:end), bOldPhi * Data(1:end-1) + bOldC, sqrt(bOldSigma2)), normpdf(Data(2:end), sOldC, sqrt(sOldSigma2))];
    case 'LN'
        sOldC = oldParam(2,2);
        sOldSigma2 = oldParam(2,3);
        Eta = [nan, nan; normpdf(Data(2:end), bOldPhi * Data(1:end-1) + bOldC, sqrt(bOldSigma2)), lognpdf(Data(2:end), sOldC, sqrt(sOldSigma2))];
    case 'P'
        %Pareto distribution
        sOldA = oldParam(2,4);
        sOldK = oldParam(2,5);
        Eta = [nan, nan; normpdf(Data(2:end), bOldPhi * Data(1:end-1) + bOldC, sqrt(bOldSigma2)), paretopdf(Data(2:end), sOldA, sOldK)];        
    case 'W'
        %Weibull distribution
        sOldA = oldParam(2,4);
        sOldB = oldParam(2,5);
        Eta = [nan, nan; normpdf(Data(2:end), bOldPhi * Data(1:end-1) + bOldC, sqrt(bOldSigma2)), weibullpdf(Data(2:end), sOldA, sOldB)];
    case 'MR'
        %mean-reverting process
        sOldPhi = oldParam(2, 1);
        sOldC = oldParam(2, 2);
        sOldSigma2 = oldParam(2, 3);        
        Eta = [nan, nan; normpdf(Data(2:end), bOldPhi * Data(1:end-1) +bOldC, sqrt(bOldSigma2)), normpdf(Data(2:end), sOldPhi * Data(1:end-1) + sOldC, sqrt(sOldSigma2))];
end
    
% Cond. prob. that observation t was generated by the mean-reverting 
% regime (first column) or the spike regime (second column) based on the 
% knowledge of data up to time t
Ksi_tt = [nan, nan; zeros(N-1, 2)];

% Cond. prob. that observation (t+1) was generated by the mean-reverting 
% regime (first column) or the spike regime (second column) based on the 
% knowledge of data up to time t
Ksi_t1t = [oldKsi_t1t_10; zeros(N-1, 2)];

% Cond. prob. that observation t was generated by the mean-reverting regime 
% (first column) or the spike regime (second column) based on the knowledge 
% of the whole data set (data up to time T), i.e. smoothed inferences
Ksi_tT = [nan, nan; zeros(N-1, 2)];

for row = 2:N
    Ksi_tt(row, :) = Eta(row, :) .* Ksi_t1t(row-1, :) / sum( Eta(row, :) .* Ksi_t1t(row-1, :) );
    LogLike(row) = log( sum( Eta(row, :) .* Ksi_t1t(row-1, :) ) );
    Ksi_t1t(row, :) = Ksi_tt(row, :) * oldP;
end

Ksi_tT(end, :) = Ksi_tt(end, :);

for row = (N-1):-1:2
    Ksi_tT(row, :) = Ksi_tt(row, :) .* ( ( Ksi_tT(row+1, :) ./ Ksi_t1t(row, :) ) * oldP' );
end
    
%Log-likelihood
Aux = Eta(2:end,:) .* Ksi_tT(2:end,:);
LogL = sum(log(Aux(:,1) + Aux(:,2)));  
    
%Transition probabilities
p_11 = sum( oldP(1,1) * Ksi_tT(3:end, 1) .* Ksi_tt(2:end-1, 1) ./ Ksi_t1t(2:end-1, 1) ) / sum( Ksi_tT(2:end-1, 1) );
p_22 = sum( oldP(2,2) * Ksi_tT(3:end, 2) .* Ksi_tt(2:end-1, 2) ./ Ksi_t1t(2:end-1, 2) ) / sum( Ksi_tT(2:end-1, 2) );
newP = [p_11, 1-p_11; 1-p_22, p_22];

newKsi_t1t_10 =  Ksi_tT(2,:);

%==========================================================================
% mrs_EM_Smoother auxiliary functions 
%==========================================================================

function Y = paretopdf(X, A, K)
if (A>0 & K>0 & sum(X<K)==0)
    Y = (A*(K^A))./(X.^(A+1));
else
    Y = nan;
end

function Y = weibullpdf(X, A, B)
if (A>0 & B>0 & sum(X<=0)==0)
    Y = A*B*(X.^(B-1)).*exp(-A*(X.^B));
else
    Y = nan;
end