/*
** Hamilton Model, mean deviation with lags
** Assuming 2-state 1st-order Markov Switching
*/
use gpe2;

@ Hamilton's Data: 1951:I to 1984:IV @
load data[]=course\gnp82.txt;
/*
@ Kim's Data @
load data[]=course\gdp4795.txt;
data=data[20:152]; @ 1951:I to 1984:IV @
*/
x=packr(100*ln(data./lagn(data,1)));
/*
compute global variable _cases for the model
_cases is used in mk1_pdf, mk1_prob, mk1_probs, mk1_llf
1st order Markov process: states=2, lags=4
*/
_cases=mk1_cases(2,4); 

call reset;
_nlopt=1;
_method=5;
_iter=200;
_print=1;
/*
_b={-0.21320,1.12828,0.08983,-0.01861,-0.17434,-0.08392,0.79619,2.20657,-1.15590};
_b={1.152,-0.3577,0.014,-0.058,-0.247,-0.213,2.0,1,2};
_b={-0.2,1.0,0.1,0.0,-0.2,-0.1,1.0,2.0,-1.0};
the list of estimated parameters is arranged in the order 
starting with the regression parameters (defined in mk1_pdf)
followed by the transition probability parameters (defined in mk1_prob) 
*/
_b={1.0,0.0,0.1,0.0,0.0,0.0,1.0,1.5,1.5}; @ 4-lag model @

call estimate(&mk1_llf,x);

call mk1_probs(x,__b);

end;

/*
Hamilton Model
based on 2-state 1st-order Markov Switching
1. Mean-deviation model
mu1=b[1], mu2=b[2]
2. with 4 autoregressive lags:
rho1=b[3], rho2=b[4], rho3=b[5], rho4=b[6]
3. constant variance:
se=b[7]
4. Transition probabilities:
prob(S(t)=1|S(t-1)=1)=b[8], prob(S(t)=2|S(t-1)=1)=1-b[8]
prob(S(t)=1|S(t-1)=2)=b[9], prob(S(t)=2|S(t-1)=2)=1-b[9]
prob={p11 p12,   where p11=b[8],   p12=b[9]
      p21 p22}         p21=1-b[8], p22=1-b[9]
*/

/* 
** user defined conditional pdf for every case
** input: x (data matrix), b (regression parameters), and cases
** output: pdf for every case and observation
** need a global variable _case:
** number of cases is states^(lags+1)
*/
proc mk1_pdf(x,b);
    local x1,x2,xx,xs,f,e,n,m;
    local bm,bs,ba,bb,r,c,i,j,k,i1;
    
    n=rows(x);         @ number of observations @
    k=rows(b);         @ number of regression parameters @
    r=rows(_cases);    @ number of cases=states^(lags+1) @
    c=cols(_cases);    @ numbers of lags+1 @
    m=int(exp(ln(r)/c));    @ number of states @
    
    bm=b[1:m];         @ m mean parameters @
    ba=b[m+1:m+c-1];   @ c-1 auto. lag parameters @
    bs=b[m+c:m+c];     @ 1 s.e. parameter @
    bb=1|-ba;          @ check: m+c=k @
    
    xx=zeros(n,m);
    j=1;
    do until j>m;
      xx[.,j]=x-b[j];
      j=j+1;
    endo;

    e=zeros(n,r);    
    i=1;
    do until i>r; @ r=states^(lags+1) @
        i1=1;
        do until i1>m; 
            if _cases[i,1]==i1; 
                xs=xx[.,i1]; 
                break;
            endif;
            i1=i1+1;
        endo;
        j=2;
        do until j>c; @ c=lags+1 @
            i1=1;
            do until i1>m; 
                if _cases[i,j]==i1; 
                    xs=xs~lagn(xx[.,i1],j-1); 
                    break;
                endif;
                i1=i1+1;
            endo;
            j=j+1;
        endo;
        e[.,i]=xs*bb;
        i=i+1;
    endo;
    bs=vec(reshape(bs,r/m,m));
    f=pdfn(packr(e)./bs')./bs';
    
    retp(f);
endp;

/* 
** define and compute probability for every case
** based on transition probabilities
** input: b (probability parameters), and cases
** output: pr case probabilities over lags+1
** need a global variable _case:
** number of cases is states^(lags+1)
*/
proc mk1_prob(b);
    local m,p,r,c,i,j,k;
    local case1,prob,a,pr,pr1,pr0;
    
    r=rows(_cases);   @ number of cases = states^(lags+1) @
    c=cols(_cases);   @ number of lags+1 @
    m=int(exp(ln(r)/c));  @ number of states @
        
    /* transition probabilities and transformation */
    p=b;   @ m*(m-1) unknown probability parameters @
    
    p=exp(b)./(1+exp(b));
    
    /*
    p=cdfn(b);
    */
        
    @ probability matrix @
    prob=reshape(p,m-1,m);
    prob=prob|(1-sumc(prob)');
    
    @ steady-state probabilities: @
    a=(eye(m)-prob)|ones(1,m);
    pr0=invpd(a'a)*a'*(zeros(m,1)|1);  @ mx1 @
    
    pr=zeros(r,c);
    i=1;
    do until i>r;
        case1=_cases[i,.]';
        pr[i,c]=pr0[case1[c]];
        j=c;
        do until j<2;
            pr[i,j-1]=prob[case1[j-1],case1[j]];
            j=j-1;
        endo;
        i=i+1;
    endo;
    
    retp(pr);
endp;

/*** do not change below ***/

proc mk1_llf(x,b);
    local m,k,k1,bp,bb,ll,lik,prob,prob0,prob1,pdfs;
    @ _cases is obtained as a global variable @
    m=int(exp(ln(rows(_cases))/cols(_cases))); 
                      @ number of states @    
    k=rows(b);        @ number of total parameters @
    k1=m*(m-1);       @ number of prabability parameters @
    bp=b[k-k1+1:k];   @ prabability parameters @
    bb=b[1:k-k1];     @ regression parameters @
    
    prob=mk1_prob(bp);   @ case probabilities @
    pdfs=mk1_pdf(x,bb);  @ case conditional pdf @
    {prob0,prob1}=mk1_probf1(pdfs,prob);
    lik=sumc((prob0.*pdfs)');
    ll=sumc(ln(lik));
    retp(ll);
endp;

proc (0) = mk1_probs(x,b);
    local m,k,k1,bp,bb,ll,lik,prob,prob0,prob1,prob2,pdfs;
    @ _cases is obtained as a global variable @
    m=int(exp(ln(rows(_cases))/cols(_cases))); 
                      @ number of states @    
    k=rows(b);        @ number of total parameters @
    k1=m*(m-1);       @ number of prabability parameters @
    bp=b[k-k1+1:k];   @ prabability parameters @
    bb=b[1:k-k1];     @ regression parameters @
    
    prob=mk1_prob(bp);   @ case probabilities @
    pdfs=mk1_pdf(x,bb);  @ case conditional pdf @
    {prob0,prob1}=mk1_probf1(pdfs,prob);
    @ lik=sumc((prob0.*pdfs)'); @
    prob2=mk1_probf2(pdfs,prob,prob0,prob1);
    
    call mk1_print(prob1,m,0); @ filtered prob @
    call mk1_print(prob2,m,0); @ smoothed prob @
endp;

/* 
** define cases based on 1-st order Markov process
** and the number of autoregressive lags
** number of cases is states^(lags+1)
*/
proc mk1_cases(m,lag);
    local i,j,c,c1,c2,cc1,cc2;
    @ m-state 1-st order Markov process @
    @ with lag autoregressive lags @

    c=seqa(1,1,m);
    c1=c'; c2=c;
    
    i=1;
    do until i>lag;
        cc1=c1;
        cc2=c2;
        j=2;
        do until j>m;        
            c1=c1|cc1;
            c2=c2|cc2;            
            j=j+1;
        endo;
        c2=vec(c1)~c2;        
        i=i+1;
    endo;
    
    retp(c2);    
endp;

/*
** compute filtered probabilities
** input: pdfs, prob
** output: prob0, prob1 or Pr(St|H(t-1)), Pr(S(t)|H(t))
*/
proc (2) = mk1_probf1(pdfs,prob);
    local pr0,n,m,c,r,i,j;
    local lik,prob0,prob1;
        
    r=rows(prob);    @ number of cases=states^(lags+1) @
    c=cols(prob);    @ number of lags+1 @
    m=int(exp(ln(r)/c));  @ number of states @
    
    c=cols(pdfs);    @ cases = rows(prob) @
    n=rows(pdfs);
    
    prob0=zeros(n,c);
    prob1=zeros(n,c);
    lik=zeros(n,1);    
    pr0=prodc(prob[.,2:cols(prob)]');
    @ product of probabilities of all previous states @
    
    i=1;
    do until i>n;
        pr0=prob[.,1].*pr0;
        prob0[i,.]=pr0';
                                
        lik[i]=pdfs[i,.]*pr0;
        pr0=(pdfs[i,.]'.*pr0)./lik[i];
        prob1[i,.]=pr0';
                
        pr0=sumc(reshape(pr0,c/m,m)');
        pr0=reshape(pr0,c,1);       
        i=i+1;
    endo;

    retp(prob0,prob1);
endp;

/*
** compute smoothed probabilities
** input: pdfs,prob,prob1
** output: prob2 or Pr(S(t)|H(N))
*/
proc mk1_probf2(pdfs,prob,prob0,prob1);    
    local i,j,k,m,n,r,c,pr1,pr2,pdfs1,prob2,lik;
    
    r=rows(prob);   @ number of cases=states^(lags+1) @
    c=cols(prob);   @ number of lags+1 @
    m=int(exp(ln(r)/c)); @ number of states @
    
    n=rows(prob1);  @ number of observations @
    c=cols(prob1);  @ cases = rows(prob) @
    lik=sumc((prob0.*pdfs)');
    prob2=prob1;    @ initialization @
    pdfs1=zeros(n,c);
    
    i=1;
    do until i>n;
        pr2=prob1[i,.]';
        
        k=1;
        do until k>c;
            pr1=zeros(c,1);
            pr1[k]=1;
        
            j=i+1;
            do until j>n;
                pr1=sumc(reshape(pr1,c/m,m)');
                pr1=reshape(pr1,c,1);
            
                pr1=prob[.,1].*pr1;            
                pdfs1[j,k]=pdfs[j,.]*pr1;
                pr1=(pdfs[j,.]'.*pr1)./pdfs1[j,k];
            
                j=j+1;       
            endo;
        
            k=k+1;
        endo;
            
        j=i+1;
        do until j>n;                   
            pr2=pr2.*(pdfs1[j,.]'./lik[j]);
            j=j+1;
        endo;
                
        prob2[i,.]=pr2';
        
        i=i+1;
    endo;
    
    retp(prob2);
endp;

/*
** compute and print current (lag=0) or lagged (lag>0)
** filtered or smoothed probabilities (pr)
*/
proc mk1_print(pr,m,lag);
    local i,j,k,n,c,pr1,prob;

    n=rows(pr);   @ number of observations @
    c=cols(pr);   @ number of cases @
                  @ ln(c)/ln(m) > lag >= 0 @
    prob=zeros(n,m);
    k=lag+1;
    i=1;
    do until i>n;
        pr1=pr[i,.]';    
        j=c/m;
        do until j<m^k;
            pr1=sumc(reshape(pr1,j,m)');
            j=j/m;
        endo;                
        do until j<m;
            pr1=sumc(reshape(pr1,m,j));
            j=j/m;
        endo;
        
        print pr1';
        prob[i,.]=pr1';
        i=i+1;
    endo;

    library pgraph;
    graphset;
    i=1;
    do until i>m;
    call xy(seqa(1,1,n),prob[.,i]);
    i=i+1;
    endo;
    retp(prob);
endp;