function xxx = lect01regest
clear;
format compact;
data = load('c:/dowjones.dat')';
N = size(data,2);
% summary stats
mm = mean(data)
ss = sqrt(mean((data-mm).^2))
%plot(1:N , data)
% initial parameters
par = [.0, .0, 0.7*ss, 1.5*ss, 0.02, 0.02];
% optimization options
optopt = optimset('MaxIter',100,'Display','iter');
[x, fval, exitflag, output, jacob, hess] = fminunc(@reglik, par, optopt, data);
par = x
stderrs = diag(sqrt(inv(hess)))'
lik =-fval
[lik,ksi] = reglikksi(par, data);
subplot(2,1,1)
plot(1:N , ksi(:,2))
subplot(2,1,2)
plot(1:N, data)


function lik = reglik(par, x)
% the likelihood based on the parameters and the data
%
% allocate parameters
m1  = par(1);
m2  = par(2);
s1  = par(3);
s2  = par(4);
p12 = par(5);
p21 = par(6);
% create densities, transition matrix
dens = [ npdf(x, m1-.5*s1^2, s1) 
         npdf(x, m2-.5*s2^2, s2) ];
cp   = [ 1-p12   p21 
          p12   1-p21  ];
x0   = [ p21
         p12 ]/(p21+p12);
% compute likelihood
[lik, ksi1] = lik(x0, dens, cp);
% to MINIMIZE
lik = -lik;

function [lik, ksi1] = reglikksi(par, x)
% the likelihood based on the parameters and the data
% output includes the filtered probabilities
% allocate parameters
m1  = par(1);
m2  = par(2);
s1  = par(3);
s2  = par(4);
p12 = par(5);
p21 = par(6);
% create densities, transition matrix
dens = [ npdf(x, m1-.5*s1^2, s1) 
         npdf(x, m2-.5*s2^2, s2) ];
cp   = [ 1-p12   p21 
          p12   1-p21  ];
x0   = [ p21
         p12 ]/(p21+p12);
% compute likelihood
[lik, ksi1] = lik(x0, dens, cp);

function y = npdf(x, m, s)
% the normal probability distribution function
y = exp(-.5 * (x-m).^2 / s^2) / s/sqrt(2*pi);

function [lik, ksi1] = lik(x0, dens, cp)
% computes the likelihood function for a regime switching model
% based on the initial distribution and the density values
%
n = size(dens,2);
m = size(dens,1);
if x0==0 % if x0=0 guess ergodic initial distribution
    cp0 = cp^64;
    x0 = cp(fix(size(cp,1)/2),:)';
end
ksi  = zeros(m,n+1); % ksi  is the matrix of the forecasted probabilities
ksi1 = zeros(m,n);   % ksi1 is the matrix of the filtered probabilities
ksi(:,1) = x0;
for indx=1:n
    zeta = ksi(:,indx) .* dens(:,indx); % weight densities
    zeta = zeta/sum(zeta);
    ksi1(:,indx) = zeta ;    
    zeta1 = cp' * zeta;
    ksi(:,indx+1) = zeta1; % update forecast
end
ksi(:,n+1) = []; % discard last forecast
lik = sum(log(abs(sum(ksi .* dens))));
ksi1 = ksi1'; % Transpose to get decent plots