📄 mrs_est.m
字号:
function [Ksi_tT,Param,P,Ksi_t1t_10,LogL]=mrs_est(Data,Model,startParam,startP,startKsi_t1t_10)
%MRS_EST Estimate 2-state (Markov) regime switching model.
% [KSI_TT,PARAM,P]=MRS_EST(DATA,MODEL) returns smoothed inferences
% KSI_TT, estimated parameters PARAM and transition matrix P of a 2-state
% (Markov) regime switching model fitted to time series DATA. The base
% regime is driven by a mean-reverting process while the spike regime by
% a Gaussian (MODEL='G'), lognormal ('LN'), Pareto ('P') or Weibull ('W')
% random variable or another mean-reverting process (MODEL='MR'). The first
% column (KSI_TT) or row (PARAM, P) contains results for the base regime
% and the second column/row for the spike regime.
% [KSI_TT,PARAM,P,KSI_T1T_10,LOGL]=MRS_EST(DATA,MODEL) additionally
% returns probabilities KSI_T1T_10 classifying the first observation to
% one of the regimes and log-likelihood LOGL of the fitted model.
% MRS_EST(DATA,MODEL,STARTPARAM,STARTP,STARTKSIT1T_10) allows to
% specify initial (model-dependent) parameter estimates STARTPARAM
% (a 1x10 vector), initial transition matrix STARTP and initial estimates
% STARTKSIT1T_10 of probabilities classifying the first observation.
% Default values for the latter three parameters are:
% switch MODEL
% case 'G'
% STARTPARAM = [.5, 3, 1, nan, nan; nan, 3, .5, nan, nan];
% case 'LN'
% STARTPARAM = [.5, 3, 1, nan, nan; nan, 3, 1, nan, nan];
% case 'P'
% STARTPARAM = [.5, 3, 1, nan, nan; nan, nan, nan, 7, .5];
% case 'W'
% STARTPARAM = [.5, 3, 1, nan, nan; nan, nan, nan, 1, 1];
% case 'MR'
% STARTPARAM = [.5, 3, 1, nan, nan; .5, 3, 1, nan, nan];
% end
% STARTP = [.5, .5; .2, .8];
% STARTKSIT1T_10 = [.6, .4];
%
% See also MRS_PLOT, MRS_SIM.
%
% Reference(s):
% [1] R.Weron (2007) 'Modeling and Forecasting Electricity Loads and
% Prices: A Statistical Approach', Wiley, Chichester.
% Written by Jakub Jurdziak and Rafal Weron (2006.09.21)
% Copyright (c) 2006 by Rafal Weron
% Set default initial estimates
if nargin < 5
startKsi_t1t_10 = [.6, .4];
end
if nargin < 4
startP = [.5, .5; .2, .8];
end
if nargin < 3
switch Model
case 'G'
%Gaussian distribution
startParam = [.5, 3, 1, nan, nan; nan, 3, .5, nan, nan];
case 'LN'
%log-normal distribution
startParam = [.5, 3, 1, nan, nan; nan, 3, 1, nan, nan];
case 'P'
%Pareto distribution
startParam = [.5, 3, 1, nan, nan; nan, nan, nan, 7, .5];
case 'W'
%Weibull distribution
startParam = [.5, 3, 1, nan, nan; nan, nan, nan, 1, 1];
case 'MR'
%mean-reverting process
startParam = [.5, 3, 1, nan, nan; .5, 3, 1, nan, nan];
end
end
% Initial estimation step
[newKsi_tT, newP, newKsi_t1t_10, newLogL] = mrs_EM_Smoother(Data, Model, startParam, startP, startKsi_t1t_10);
newParam = mrs_EM_MLE(Data, Model, newKsi_tT, startParam);
% Estimation loop
iteration = 1;
diff = 1;
while and(diff > 10^(-8), iteration < 100)
oldP = newP;
oldKsi_t1t_10 = newKsi_t1t_10;
oldParam = newParam;
oldKsi_tT = newKsi_tT;
oldLogL = newLogL;
[newKsi_tT, newP, newKsi_t1t_10, newLogL] = mrs_EM_Smoother(Data, Model, oldParam, oldP, oldKsi_t1t_10);
newParam = mrs_EM_MLE(Data, Model, newKsi_tT, oldParam);
diffP = max(abs(newP-oldP));
diffKsi_t1t_10 = max(abs(newKsi_t1t_10-oldKsi_t1t_10));
diffParam = max(abs(newParam-oldParam));
diffKsi_tT = max(abs(newKsi_tT-oldKsi_tT));
diffLogL = abs(newLogL-oldLogL);
diff = max([diffP diffKsi_t1t_10 diffParam diffKsi_tT diffLogL]);
iteration = iteration+1;
end
% Results
Ksi_tT = newKsi_tT;
Param = newParam;
P = newP;
Ksi_t1t_10 = newKsi_t1t_10;
LogL = newLogL;
Summary = mrs_Summary(Model, Param, P, LogL);
switch Model
case 'G'
% Gaussian distribution
disp([' ']);
disp(['Two state regime switching model with Gaussian spike distribution']);
disp([' ']);
disp(['regime beta_i c_i sigma^2_i E(Y_{t,i}) Var(Y_{t,i}) q_ii P(R=i)']);
disp(['base ' sprintf('%7.5f',Summary(1,1)) ' ' sprintf('%7.5f',Summary(1,2)) ' ' sprintf('%7.5f',Summary(1,3)) ' ' ...
sprintf('%7.5f',Summary(1,4)) ' ' sprintf('%7.5f',Summary(1,5)) ' ' sprintf('%7.5f',Summary(1,6)) ' ' sprintf('%7.5f',Summary(1,7))]);
disp(['spike - ' sprintf('%7.5f',Summary(2,2)) ' ' sprintf('%7.5f',Summary(2,3)) ' ' ...
sprintf('%7.5f',Summary(2,4)) ' ' sprintf('%7.5f',Summary(2,5)) ' ' sprintf('%7.5f',Summary(2,6)) ' ' sprintf('%7.5f',Summary(2,7))]);
case 'LN'
% log-normal distribution
disp([' ']);
disp(['Two state regime switching model with lognormal spike distribution']);
disp([' ']);
disp(['regime beta_i c_i sigma^2_i E(Y_{t,i}) Var(Y_{t,i}) q_ii P(R=i)']);
disp(['base ' sprintf('%7.5f',Summary(1,1)) ' ' sprintf('%7.5f',Summary(1,2)) ' ' sprintf('%7.5f',Summary(1,3)) ' ' ...
sprintf('%7.5f',Summary(1,4)) ' ' sprintf('%7.5f',Summary(1,5)) ' ' sprintf('%7.5f',Summary(1,6)) ' ' sprintf('%7.5f',Summary(1,7))]);
disp(['spike - ' sprintf('%7.5f',Summary(2,2)) ' ' sprintf('%7.5f',Summary(2,3)) ' ' ...
sprintf('%7.5f',Summary(2,4)) ' ' sprintf('%7.5f',Summary(2,5)) ' ' sprintf('%7.5f',Summary(2,6)) ' ' sprintf('%7.5f',Summary(2,7))]);
case 'P'
% Pareto distribution
disp([' ']);
disp(['Two state regime switching model with Pareto ditributed spikes']);
disp([' ']);
disp(['regime beta_i c_i / a s^2_i / k E(Y_{t,i}) Var(Y_{t,i}) q_ii P(R=i)']);
disp(['base ' sprintf('%7.5f',Summary(1,1)) ' ' sprintf('%7.5f',Summary(1,2)) ' ' sprintf('%7.5f',Summary(1,3)) ' ' ...
sprintf('%7.5f',Summary(1,4)) ' ' sprintf('%7.5f',Summary(1,5)) ' ' sprintf('%7.5f',Summary(1,6)) ' ' sprintf('%7.5f',Summary(1,7))]);
disp(['spike - ' sprintf('%7.5f',Summary(2,2)) ' ' sprintf('%7.5f',Summary(2,3)) ' ' ...
sprintf('%7.5f',Summary(2,4)) ' ' sprintf('%7.5f',Summary(2,5)) ' ' sprintf('%7.5f',Summary(2,6)) ' ' sprintf('%7.5f',Summary(2,7))]);
case 'W'
% Weibull distribution
disp([' ']);
disp(['Two state regime switching model with Weibull ditributed spikes']);
disp([' ']);
disp(['regime beta_i c_i / a s^2_i / b E(Y_{t,i}) Var(Y_{t,i}) q_ii P(R=i)']);
disp(['base ' sprintf('%7.5f',Summary(1,1)) ' ' sprintf('%7.5f',Summary(1,2)) ' ' sprintf('%7.5f',Summary(1,3)) ' ' ...
sprintf('%7.5f',Summary(1,4)) ' ' sprintf('%7.5f',Summary(1,5)) ' ' sprintf('%7.5f',Summary(1,6)) ' ' sprintf('%7.5f',Summary(1,7))]);
disp(['spike - ' sprintf('%7.5f',Summary(2,2)) ' ' sprintf('%7.5f',Summary(2,3)) ' ' ...
sprintf('%7.5f',Summary(2,4)) ' ' sprintf('%7.5f',Summary(2,5)) ' ' sprintf('%7.5f',Summary(2,6)) ' ' sprintf('%7.5f',Summary(2,7))]);
case 'MR'
% mean-reverting process
disp([' ']);
disp(['Two state regime switching model with mean-reverting process for spikes']);
disp([' ']);
disp(['regime beta_i c_i sigma^2_i E(Y_{t,i}) Var(Y_{t,i}) q_ii P(R=i)']);
disp(['base ' sprintf('%7.5f',Summary(1,1)) ' ' sprintf('%7.5f',Summary(1,2)) ' ' sprintf('%7.5f',Summary(1,3)) ' ' ...
sprintf('%7.5f',Summary(1,4)) ' ' sprintf('%7.5f',Summary(1,5)) ' ' sprintf('%7.5f',Summary(1,6)) ' ' sprintf('%7.5f',Summary(1,7))]);
disp(['spike ' sprintf('%7.5f',Summary(2,1)) ' ' sprintf('%7.5f',Summary(2,2)) ' ' sprintf('%7.5f',Summary(2,3)) ' ' ...
sprintf('%7.5f',Summary(2,4)) ' ' sprintf('%7.5f',Summary(2,5)) ' ' sprintf('%7.5f',Summary(2,6)) ' ' sprintf('%7.5f',Summary(2,7))]);
end
%=========================================================================
% Internally used routine(s)
%=========================================================================
%=========================================================================
% mrs_Summary
%=========================================================================
function Summary = mrs_Summary(Model, Param, P, logL)
%mean and variance for base regime
bParam = Param(1,:);
bPhi = bParam(1);
bC = bParam(2);
bSigma2 = bParam(3);
bMean = bC / (1-bPhi);
bVar = bSigma2 / (1-bPhi^2);
%mean and variance for spike regime
sParam = Param(2,:);
[sMean, sVar] = mrs_Moments(Model, sParam);
%unconditional probabilities
p_11 = P(1,1);
p_22 = P(2,2);
P_Rt_1 = (1-p_22) / (2-p_11-p_22);
P_Rt_2 = (1-p_11) / (2-p_11-p_22);
%summary matrix
switch Model
case 'G'
%Gaussian distribution
sC = sParam(2);
sSigma2 = sParam(3);
Summary = [ 1-bPhi, bC, bSigma2, bMean, bVar, p_11, P_Rt_1;
nan, sC, sSigma2, sMean, sVar, p_22, P_Rt_2;
logL, nan, nan, nan, nan, nan, nan];
case 'LN'
sC = sParam(2);
sSigma2 = sParam(3);
Summary = [ 1-bPhi, bC, bSigma2, bMean, bVar, p_11, P_Rt_1;
nan, sC, sSigma2, sMean, sVar, p_22, P_Rt_2;
logL, nan, nan, nan, nan, nan, nan];
case 'P'
%Pareto distribution
sA = sParam(4);
sK = sParam(5);
Summary = [ 1-bPhi, bC, bSigma2, bMean, bVar, p_11, P_Rt_1;
nan, sA, sK, sMean, sVar, p_22, P_Rt_2;
logL, nan, nan, nan, nan, nan, nan];
case 'W'
%Weibull distribution
sA = sParam(4);
sB = sParam(5);
Summary = [ 1-bPhi, bC, bSigma2, bMean, bVar, p_11, P_Rt_1;
nan, sA, sB, sMean, sVar, p_22, P_Rt_2;
logL, nan, nan, nan, nan, nan, nan];
case 'MR'
%mean-reverting process
sPhi = sParam(1);
sC = sParam(2);
sSigma2 = sParam(3);
Summary = [ 1-bPhi, bC, bSigma2, bMean, bVar, p_11, P_Rt_1;
1-sPhi, sC, sSigma2, sMean, sVar, p_22, P_Rt_2;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -