/*
** Mixture of Two Normal Distributions
** Stochastic switching model, based on Lesson 6.5
** See Greene [1999], Chapter 4
*/
use gpe2;

load data[21,2]=gpe\yed20.txt;
x=data[2:21,1]/10; @ income data: scaling may help @
b1={3,2};
b2={3,2};
b3={3,2};

call reset;
_nlopt=2;
_method=5;
_iter=100;

@ 2-state model @
_b=b1|b2|0;
call estimate(&llf,x);
print prob(x,__b);
/*
@ 3-state model @
_b=b1|b2|b3|0|0;
call estimate(&llf,x);
print prob(x,__b);
*/
end;

/* 
** user needs to define pdf for all states
** number of states = rows of returned pdf matrix
** 2-state model:
mixture of two normal distributions
mu1=b[1], mu2=b[3]
se1=b[2], se2=b[4]
prob.(drawn from the 1st distribution)=b[5]
** 3-state model:
mixture of three normal distributions
mu1=b[1], mu2=b[3], mu3=b[5]
se1=b[2], se2=b[4], se2=b[6]
prob.(drawn from the 1st distribution)=b[7]
prob.(drawn from the 2nd distribution)=b[8]
*/

proc pdf(x,b); @ 2-state pdfs @
    local pdf1,pdf2;
    pdf1=pdfn((x-b[1])/b[2])/b[2];
    pdf2=pdfn((x-b[3])/b[4])/b[4];
    retp(pdf1~pdf2);
endp;
/*
proc pdf(x,b);  @ 3-state pdfs @
    local pdf1,pdf2,pdf3;
    pdf1=pdfn((x-b[1])/b[2])/b[2];
    pdf2=pdfn((x-b[3])/b[4])/b[4];
    pdf3=pdfn((x-b[5])/b[6])/b[6];
    retp(pdf1~pdf2~pdf3);
endp;
*/
/*
** component log-likelihood function
*/
proc llf(x,b);
    local pdfs,ll,p,m,k;
    pdfs=pdf(x,b);
    m=cols(pdfs);       @ number of states @
    k=rows(b);          @ number of switches = m-1 @
    /*
    p=b[k-m+2:k];
    */
    p=cdfn(b[k-m+2:k]);
    p=p|(1-sumc(p));
    ll=pdfs*p;
    ll=ln(ll.*(ll.>0)+(0.0)*(ll.<=0));
    retp(ll);
endp;

/*
** conditional probability of each state
*/
proc prob(x,b);
    local pdfs,p,m,k;
    pdfs=pdf(x,b);      @ get pdf for all states @
    m=cols(pdfs);       @ number of states @
    k=rows(b);
    /*
    p=b[k-m+2:k];
    */
    p=cdfn(b[k-m+2:k]);
    p=p|(1-sumc(p));
    retp((p'.*pdfs)./(pdfs*p));
endp;