function SimVal=mrs_sim(Model,Param,Ksi_tT,Data1,simNo,Stats)
%MRS_SIM Simulate trajectories of a 2-state (Markov) regime switching model. 
%   MRS_SIM(MODEL,PARAM,KSI_TT,DATA1,SIMNO) generates SIMNO trajectories of  
%   the 2-state (Markov) regime switching model with Gaussian (MODEL='G'),
%   lognormal ('LN'), Pareto ('P') or Weibull ('W') distributed spikes or 
%   a mean-reverting process for the spike regime (MODEL='MR'). PARAM and
%   KSI_TT are the model parameters and smoothed inferences, respectively, 
%   returned by MRS_EST. DATA1 is the starting point for the simulated
%   trajectories.
%   MRS_SIM(MODEL,PARAM,KSI_TT,DATA1,SIMNO,1) additionally displays  
%   statistics for the first differences of the simulated trajectories.
%   SIMVAL=MRS_SIM(MODEL,PARAM,KSI_TT,DATA1,SIMNO) returns a matrix of 
%   simulated trajectories.
%
%   See also MRS_EST, MRS_PLOT.
%
%   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


if nargin<6
    Stats = 0;
end;

% Define length of the simulated trajectories
simLen = length(Ksi_tT)-1;

% Discard the first (NaN) smoothed inference
Ksi_tT = Ksi_tT(2:end,1);

% Generate a matrix of random numbers for specifying regimes 
rndRegime = rand(simLen,simNo);
bRegime = rndRegime < repmat(Ksi_tT, 1, simNo);
sRegime = 1-bRegime;

% Spike regime simulation
switch Model
    case 'G'
        sSimVal = normrnd(Param(2,2), sqrt(Param(2,3)), simLen, simNo);      
    case 'LN'
        sSimVal = lognrnd(Param(2,2), sqrt(Param(2,3)), simLen, simNo);
    case 'P'    
        sSimVal = paretornd(Param(2,4), Param(2,5), simLen, simNo);
    case 'W'    
        sSimVal = weibullrnd(Param(2,4), Param(2,5), simLen, simNo);
    case 'MR' 
        sSimVal = Param(2,1)*Data1 + normrnd(Param(2,2), sqrt(Param(2,3)) ,1 , simNo);
        for i=2:simLen
            sSimVal(i,:) = Param(2,1) * sSimVal(i-1,:) + normrnd(Param(2,2), sqrt(Param(2,3)) ,1 , simNo);
        end
end

% Base regime simulation
bSimVal = Param(1,1)*Data1 + normrnd(Param(1,2), sqrt(Param(1,3)) ,1 , simNo);
for i=2:simLen
    bSimVal(i,:) = Param(1,1) * bSimVal(i-1,:) + normrnd(Param(1,2), sqrt(Param(1,3)) ,1 , simNo);
end  

SimVal = bSimVal .* bRegime + sSimVal .* sRegime;

% Statistics for the first differences of the simulated trajectories
if Stats,
    DSimVal = diff(SimVal);
    meanS = mean(DSimVal);
    varS = var(DSimVal);
    meanvS = repmat(meanS,length(Ksi_tT)-1,1);
    skewS = mean((DSimVal - meanvS) .^ 3 ) ./ (varS.^ 1.5);
    kurtS = -3 + mean((DSimVal - meanvS) .^ 4) ./ (varS.^ 2);
    minS = min(DSimVal);
    maxS = max(DSimVal);
    uJumpsS = sum(DSimVal>0.3);
    dJumpsS = sum(DSimVal<-0.3);
    uEJumpsS = sum((DSimVal>0.3).*DSimVal) ./ uJumpsS;
    dEJumpsS = sum((DSimVal<-0.3).*DSimVal) ./ dJumpsS;
    q99S = prctile(DSimVal,99);
    q995S = prctile(DSimVal,99.5);

    disp([' ']);        
    disp(['Mean values of descriptive statistics for differenced trajectories']);                        
    switch Model
        case 'G'
            %Gaussian distribution
            disp(['  - 2-state regime switching model with Gaussian spikes:']);
        case 'LN'
            %log-normal distribution
            disp(['  - 2-state regime switching model with lognormal spikes:']);        
        case 'P'
            %Pareto distribution
            disp(['  - 2-state regime switching model with Pareto spikes:']);                
        case 'W'
            %Weibull distribution
            disp(['  - 2-state regime switching model with Weibull spikes:']);                
        case 'MR'
            %mean-reverting process
            disp(['  - 2-state regime switching model with a mean-reverting process for spikes:']);
    end
    disp([' ']);
    disp(['Minimum:                 ' sprintf('%7.5f',mean(minS))]);   
    disp(['Maximum:                 ' sprintf('%7.5f',mean(maxS))]);   
    disp(['Mean:                    ' sprintf('%7.5f',mean(meanS))]);   
    disp(['Variance:                ' sprintf('%7.5f',mean(varS))]);   
    disp(['Skewness:                ' sprintf('%7.5f',mean(skewS))]);   
    disp(['Kurtosis:                ' sprintf('%7.5f',mean(kurtS))]);   
    disp(['99% quantile:            ' sprintf('%7.5f',mean(q99S))]);   
    disp(['99.5% quantile:          ' sprintf('%7.5f',mean(q995S))]);
    disp(['No. upward jumps:        ' sprintf('%7.5f',mean(uJumpsS))]);
    disp(['No. downward jumps:      ' sprintf('%7.5f',mean(dJumpsS))]);
    disp(['Mean upward jump size:   ' sprintf('%7.5f',mean(uEJumpsS))]);
    disp(['Mean downward jump size: ' sprintf('%7.5f',mean(dJumpsS))]);
end;

%=========================================================================
% Internally used routine(s)
%=========================================================================

function Y = paretornd(A, K, rowsNo, colNo)
Y = K*(1./(rand(rowsNo, colNo))).^(1/A);

function Y = weibullrnd(A, B, rowsNo, colNo)
Y = (log(rand(rowsNo,colNo))/(-A)).^(1/B);