Problem with

You are missing some code. lik_fcn is a likelihood procedure written by whoever you got the code from. It is not a built-in GAUSS function.

The code for tne lik_fcn procedure will start something like this:

proc (1) lik_fcn(x);
    local ...

You are right! I missed some code. The full code is:

new;
cls;
 
library optmum,PGRAPH;
 
format /m1 /rd 8,6;
load data[134,3]=Pib.txt; // Dados mensais: PIB e choques monetarios
// 1995.7 -- 2006.8
yy_d=(data[1:134,.]); //recebe PIB
 
t0=rows(yy_d); //numero de observacoes
 
LAG_AR=7; //numero de lags AR(p)
NO_ST=lag_ar+1; // Número de Estados
DMNSION=2^NO_ST;
 
st_mat=zeros(DMNSION,NO_ST); //Cria matriz de zeros com DMNSION linhas e NO_ST colunas
//talvez seja um passo necessario pro metodo do Kim (1994)
 
j=1;
st7=0;do until st7>1;
st6=0;do until st6>1;
st5=0;do until st5>1;
st4=0;do until st4>1;
st3=0; do until st3>1;
st2=0; do until st2>1;
st1=0; do until st1>1;
st=0; do until st>1;
 
st_mat[j,.]=st7~st6~st5~st4~st3~st2~st1~st;
 
j=j+1;
st=st+1;endo;
st1=st1+1;endo;
st2=st2+1;endo;
st3=st3+1;endo;
st4=st4+1;endo;
st5=st5+1;endo;
st6=st6+1;endo;
st7=st7+1;endo;
 
yy=yy_d[lag_ar+1:t0,1]; //Esta separando um vetor sem as 7 primeiras observacoes
 
x_mat=
yy_d[Lag_ar:T0-1,1] ~ yy_d[lag_ar-1:T0-2,1]~
yy_d[lag_ar-2:T0-3,1]~ yy_d[lag_ar-3:T0-4,1]~
yy_d[lag_ar-4:T0-5,1]~yy_d[lag_ar-5:T0-6,1]~yy_d[lag_ar-6:T0-7,1] ;
 
z_mat=
yy_d[Lag_ar:T0-1,2] ~ yy_d[lag_ar-1:T0-2,2]~
yy_d[lag_ar-2:T0-3,2]~ yy_d[lag_ar-3:T0-4,2]~
yy_d[lag_ar-4:T0-5,2]~ yy_d[lag_ar-5:T0-6,2]~
yy_d[lag_ar-6:T0-7,2]~
yy_d[Lag_ar:T0-1,3]~ yy_d[lag_ar-1:T0-2,3]~
yy_d[lag_ar-2:T0-3,3]~ yy_d[lag_ar-3:T0-4,3]~
yy_d[lag_ar-4:T0-5,3]~yy_d[lag_ar-5:T0-6,3]~
yy_d[lag_ar-6:T0-7,3];
 
Dados=yy~x_mat~z_mat;
format /m1/rd 8,4;
//"Dados sao:"; dados;
//"";
//output file= C:\Gauss\Tese\Dados.txt reset;
//dados;
output off;
 
T=rows(yy);
 
// Valores iniciais para rodar a maximização de MV
//================= Inicialização================================
 
START=2; //1995:8.....
 
PRMTR_IN={2.818921 1.035 -1.226890 0.549770 -0.237239 -0.271351 -0.306826 -0.123236
-0.230071 -0.090165 0.000000 -2.295783 -4.612134 2.965306 -0.248680 -4.764720
-0.022357 3.363684 -2.638171 -5.761053 -0.101388 3.441007 -1.235622 0.000000
0.000000 -0.570660 -0.969447 -6.304970 -0.772019 1.234718 -3.516800 -7.764870
1.826352 -3.444785 -3.323838 2.707739 1.130373 0.000000 0.000000 1.131948};
 
PRMTR_IN=PRMTR_IN';
 
// Estimação por ML
//==================================================
output off;
{xout,fout,gout,cout}=optmum(&lik_fcn,PRMTR_in);
 
PRM_FNL=TRANS(xout); // Estimativas dos parâmetros
AIC = -2*(-fout)/(T) + 2*rows(prmtr_in)/(T);
SIC = -2*(-fout)/(T) + rows(prmtr_in)*ln(T)/(T);
 
output file=c:\gauss10\Markov\pibout.txt reset;
output off;
 
//Calcular os estimadores de MV requer averiguar se estamos no ponto de máximo,
//por isso tem o cálculo do hessiano
 
"==Resultados Finais========================================================";
 
"Valor inicial dos parâmetros";
trans(prmtr_in)';
 
"Valor da Loglik ";; -1*fout;
"SIC is:";; SIC;
"AIC is:";; AIC;
 
"Parâmetros estimados:";
prm_fnl';
xout;
output off;
 
"Calculando matrix Hessiana";
hout0=hessp(&lik_fcn,xout);
hout=inv(hout0);
 
grdn_fnl=gradfd(&TRANS,xout);
Hsn_fnl=grdn_fnl*hout*grdn_fnl';
SD_fnl =sqrt(diag(Hsn_fnl)); //Standard errors of the estimated coefficients
tstudent=prm_fnl./sd_fnl;
"T-Student dos parâmetros:"; tstudent;
nobs=T;
nbparam=rows(prm_fnl);
ddl=T-nbparam;
 
"Graus de liberdade" ddl;
pvalue=2*cdftc(abs(tstudent),ddl);
 
output off;
"Erro padrão dos parâmetros:"; sd_fnl';
"T-Student dos parâmetros:"; tstudent;
"Numero de observacoes:"; nobs;
"Graus de liberdade:"; ddl;
"pvalues:";pvalue;
"===============================================================";
output off;
 
/*Salvando os parametros, desvios-padrao e p-value*/
 
param_sd=prm_fnl~sd_fnl~pvalue;
"param_sd são:"; param_sd;
"";
format /rd 8,8;
//output file= C:\Gauss\Tese\Param_sd.txt reset; // Alterar diretorio
//param_sd;
output off;
 
/*Salvando a matriz de covariancia dos parametros*/
 
Covparam = Hsn_fnl;
format /rd 8,8;
"Covparam dos parâmetros:"; covparam;
"";
//output file= C:\Gauss\Tese\Covparam.txt reset;
//covparam;
output off;
 
//Depois de determinar o melhor modelo, verificar se os choques são assimétricos
 
/*Testes de Wald*/
 
/*Wald 1: O efeito de choques negativos é nulo*/
R1={0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0};
q1={0 0};
m1=R1*prm_fnl - q1';
varm1 = R1*covparam*R1';
wald1 = m1'*inv(varm1)*m1;
 
/*Wald 2: O efeito de choques positivos é nulo*/
R2={0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0};
q2={0 0};
m2=R2*prm_fnl - q2';
varm2 = R2*covparam*R2';
wald2 = m2'*inv(varm2)*m2;
 
/*Wald 3: Teste de assimetria*/
R3={0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 -1 0 -1 0 -1 0 -1 0 -1 0 -1 0 -1 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 -1 0 -1 0 -1 0 -1 0 -1 0 -1 0 -1 0};
q3={0 0};
m3=R3*prm_fnl - q3';
varm3 = R3*covparam*R3';
wald3 = m3'*inv(varm3)*m3;
 
/*Wald 4: O efeito de choques negativos é nulo apenas no regime 1*/
R4={0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0};
q4={0};
m4=R4*prm_fnl - q4';
varm4 = R4*covparam*R4';
wald4 = m4'*inv(varm4)*m4;
 
/*Wald 5: O efeito de choques negativos é nulo apenas no regime 2*/
R5={0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0};
q5={0};
m5=R5*prm_fnl - q5';
varm5 = R5*covparam*R5';
wald5 = m5'*inv(varm5)*m5;
 
/*Wald 6: O efeito de choques positivos é nulo apenas no regime 1*/
R6={0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 0};
q6={0};
m6=R6*prm_fnl - q6';
varm6 = R6*covparam*R6';
wald6 = m6'*inv(varm6)*m6;
 
/*Wald 7: O efeito de choques positivos é nulo apenas no regime 2*/
R7={0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0};
q7={0};
m7=R7*prm_fnl - q7';
varm7 = R7*covparam*R7';
wald7 = m7'*inv(varm7)*m7;
 
/*Wald 8: Simetria entre os choques apenas no regime 1*/
R8={0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 -1 0 -1 0 -1 0 -1 0 -1 0 -1 0 -1 0 0};
q8={0};
m8=R8*prm_fnl - q8';
varm8 = R8*covparam*R8';
wald8 = m8'*inv(varm8)*m8;
 
/*Wald 9: Simetria entre os choques apenas no regime 2*/
R9={0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 -1 0 -1 0 -1 0 -1 0 -1 0 -1 0 -1 0};
q9={0};
m9=R9*prm_fnl - q9';
varm9 = R9*covparam*R9';
wald9 = m9'*inv(varm9)*m9;
 
/*Wald 10: Simetria entre choques negativos nos regimes 1 e 2*/
R10={0 0 0 0 0 0 0 0 0 0 0 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0};
q10={0};
m10=R10*prm_fnl - q10';
varm10 = R10*covparam*R10';
wald10 = m10'*inv(varm10)*m10;
 
/*Wald 11: Simetria entre choques positivos nos regimes 1 e 2*/
R11={0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 0};
q11={0};
m11=R11*prm_fnl - q11';
varm11 = R11*covparam*R11';
wald11 = m11'*inv(varm11)*m11;
 
/*Wald 12: Igualdade entre os parametros dos choques positivos e dos choques negativos nos regimes 1 e 2*/
R12={0 0 0 0 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1 0};
q12={0 0 0 0 0 0 0 0 0 0 0 0 0 0};
m12=R12*prm_fnl - q12';
varm12 = R12*covparam*R12';
wald12 = m12'*inv(varm12)*m12;
 
/*Wald 13: Igualdade entre os parametros dos choques negativos nos regimes 1 e 2 - Assimetria dos choques negativos entre regimes */
R13={0 0 0 0 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0};
q13={0 0 0 0 0 0 0};
m13=R13*prm_fnl - q13';
varm13=R13*covparam*R13';
wald13=m13'*inv(varm13)*m13;
 
/*Wald 14: Igualdade entre os parametros dos choques positivos nos regimes 1 e 2 - Assimetria dos choques positivos entre regimes */
R14={0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1 0};
q14={0 0 0 0 0 0 0};
m14=R14*prm_fnl - q14';
varm14 = R14*covparam*R14';
wald14 = m14'*inv(varm14)*m14;
 
/*Wald 15: Igualdade entre os parametros dos choques positivos e negativos nos regimes 1 e 2 */
R15={0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0};
q15={0 0 0 0 0 0 0 0 0 0 0 0 0 0};
m15=R15*prm_fnl - q15';
varm15 = R15*covparam*R15';
wald15 = m15'*inv(varm15)*m15;
 
/*Wald 16: Igualdade entre os parametros dos choques positivos e negativos no regime 1 */
R16={0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0};
q16={0 0 0 0 0 0 0};
m16=R16*prm_fnl - q16';
varm16=R16*covparam*R16';
wald16=m16'*inv(varm16)*m16;
 
/*Wald 17: Igualdade entre os parametros dos choques positivos e negativos no regime 2 */
R17={ 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0};
q17={0 0 0 0 0 0 0};
m17=R17*prm_fnl - q17';
varm17 = R17*covparam*R17';
wald17 = m17'*inv(varm17)*m17;
 
/*Wald 18: Coeficientes dos choques negativos no regime 1 são conjuntamente iguais a zero */
R18={0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0};
q18={0 0 0 0 0 0 0};
m18=R18*prm_fnl - q18';
varm18=R18*covparam*R18';
wald18=m18'*inv(varm18)*m18;
 
/*Wald 19: Coeficientes dos choques negativos no regime 2 são conjuntamente iguais a zero */
R19={0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0};
q19={0 0 0 0 0 0 0};
m19=R19*prm_fnl - q19';
varm19=R19*covparam*R19';
wald19=m19'*inv(varm19)*m19;
 
/*Wald 20: Coeficientes dos choques positivos no regime 1 são conjuntamente iguais a zero */
R20={0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0};
q20={0 0 0 0 0 0 0};
m20=R20*prm_fnl - q20';
varm20 = R20*covparam*R20';
wald20 = m20'*inv(varm20)*m20;
 
/*Wald 21: Coeficientes dos choques positivos no regime 2 são conjuntamente iguais a zero */
R21={0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0};
q21={0 0 0 0 0 0 0};
m21=R21*prm_fnl - q21';
varm21 = R21*covparam*R21';
wald21 = m21'*inv(varm21)*m21;
 
/*Wald 22: Igualdade entre os parametros dos choques positivos no regime 1 e negativos no regime 2 */
R22={0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0};
q22={0 0 0 0 0 0 0};
m22=R22*prm_fnl - q22';
varm22=R22*covparam*R22';
wald22 = m22'*inv(varm22)*m22;
 
/*Wald 23: Simetria entre choques negativos no regime 2 e choques positivos no regime 1*/
R23={0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 1 0 1 0 1 0 1 -1 0 -1 0 -1 0 -1 0 -1 0 -1 0 -1 0 0};
q23={0};
m23=R23*prm_fnl - q23';
varm23 = R23*covparam*R23';
wald23 = m23'*inv(varm23)*m23;
 
/*Wald 24: Simetria entre choques negativos no regime 1 e choques positivos no regime 2*/
R24={0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 0 -1 0 -1 0 -1 0 -1 0 -1 0 -1 0 -1 0};
q24={0};
m24=R24*prm_fnl - q24';
varm24 = R24*covparam*R24';
wald24 = m24'*inv(varm24)*m24;
 
Wald=wald1~wald2~wald3~wald4~wald5~wald6~wald7~wald8~wald9~wald10~wald11~wald12~wald13~wald14~wald15~wald16~wald17~wald18~wald19~wald20~wald21~wald22~wald23~wald24;
format /rd 8,8;
"Wald:"; wald;
"";
//output file= C:\Gauss\Tese\Wald.txt reset;
//wald;
output off;
 
/* Probabilidades Filtradas, previstas e Suavizadas*/
 
{pr_tt0,pr_tl0,out_mat}=FILTER(XOUT); //Pr[S_t=0|Y_t] and Pr[S_t=0|Y_{t-1}]
 
smooth0=SMOOTH(pr_tt0,pr_tl0); //Pr[S_t=0|Y_T], Smoothed probabilities
 
RESIDUOS=out_mat;
format /m1 /rd 8,8;
//output file=C:\Gauss\Tese\RESIDUOS.txt reset; RESIDUOS;
output off;
 
FLTR=pr_tt0~pr_tl0~smooth0;
format /m1 /rd 8,8;
//output file=C:\Gauss\Tese\Probabilidades.txt reset; FLTR;
output on;
 
/*Salvando as probabilidade filtradas*/
Filtrada=pr_tt0~1-pr_tt0;
format /m1/rd 8,8;
"filtradas are:"; filtrada;
"";
//output file= C:\Gauss\Tese\ProbFiltrada.txt reset;
//filtrada;
output off;
 
/*Salvando as probabilidade suavizadas*/
Suavizada=smooth0~1-smooth0;
format /m1/rd 8,8;
"Probabilidades Suavizadas sao:"; suavizada;
"";
//output file= C:\Gauss\Tese\ProbSuavizada.txt reset;
//suavizada;
output off;
 
/*Calculando o criterio MSC*/
Lambda=2;
K=22; /*K=((rows(prmtr_in))-3)/2;*/
T1=sumc(smooth0);
T2=sumc(1-smooth0);
MSC = -2*(-fout) + ((T1*(T1+lambda*K))/(T1-lambda*K-2)) + ((T2*(T2+lambda*K))/(T2-lambda*K-2));
Criterios=AIC~SIC~MSC;
format /m1/rd 8,8;
"Criterios de informacao sao:"; criterios;
"";
//output file= C:\Gauss\Tese\Criterios.txt reset;
//criterios;
output off;
 
xy(seqa(1,1,rows(pr_tt0)),pr_tt0);
xy(seqa(1,1,rows(smooth0)),smooth0);
xy(seqa(1,1,rows(out_mat)),out_mat~yy);
end;
 
// Fim do programa principal
 
//========================================================================
 
//========================================================================
PROC LIK_FCN(PRMTR1);
local prmtr, ppr,qpr,prob__0,prob__1,QQ, lik, j_iter, pr__0_l,pr__1_l,
F_cast, var_L,pr_vl, pr_val,likv,phi,PSIC,PSIX,
vecp,st,st1,st2,st3,st4,ST5,ST6,ST7,ST8,ST9,ST10,ST11,ST12,ST13,
pr_tr,pr_trf1,pr_trf,prob__t,prob__,pro_,j,psi1,psi2,var_c,
delta0,DELTA1,MU0,MU1, st_k,st_k1,st_k2,st_k3,st_k4,
f_cast1,f_cast2,PR_VL1,PR_VL2,pr_trf7,pr_trf0,
PR_TRF2,PR_TRF3,PR_TRF4,PR_TRF5,PR_TRF6,psiL,
TMPRY1,TMPRY2,SM_PRL,TMP_P0,SM_PR0,JJJ,MU_MAT,D_MAT,FLT_PRN,
F1,F2,TMP00,TMP0,SM_PR00,prob_dd,VAR,A,EN,BETANEG10,BETANEG11,
BETANEG20,BETANEG21,BETANEG30,BETANEG31,BETANEG40,BETANEG41,BETANEG50,BETANEG51,BETANEG60,BETANEG61,BETANEG70,BETANEG71,BETANEG80,BETANEG81,BETANEG90,BETANEG91,
BETANEG1_MAT,BETANEG2_MAT,BETANEG3_MAT,BETANEG4_MAT,BETANEG5_MAT,BETANEG6_MAT,BETANEG7_MAT,BETANEG8_MAT,BETANEG9_MAT,
BETAPOS10,BETAPOS11,BETAPOS20,BETAPOS21,BETAPOS30,BETAPOS31,BETAPOS40,BETAPOS41,BETAPOS50,BETAPOS51,BETAPOS60,BETAPOS61,
BETAPOS70,BETAPOS71,BETAPOS80,BETAPOS81,BETAPOS90,BETAPOS91,BETAPOS1_MAT,BETAPOS2_MAT,BETAPOS3_MAT,BETAPOS4_MAT,BETAPOS5_MAT,BETAPOS6_MAT,BETAPOS7_MAT,BETAPOS8_MAT,BETAPOS9_MAT;
 
PRMTR=TRANS(PRMTR1);
 
LOCATE 15,1; //PRMTR';
 
/* Definição dos parâmetros */
 
PPR=PRMTR[1,1]; //Pr[St=1/St-1=1] St=1 indica expansao
QPR=PRMTR[2,1]; //Pr[St=0/St-1=0] St=0 indica recessao
 
MU0=PRMTR[3,1];
MU1=PRMTR[4,1];
PHI=PRMTR[5:11,1];
 
BETANEG10=PRMTR[12,1];
BETANEG11=PRMTR[13,1];
BETANEG20=PRMTR[14,1];
BETANEG21=PRMTR[15,1];
BETANEG30=PRMTR[16,1];
BETANEG31=PRMTR[17,1];
BETANEG40=PRMTR[18,1];
BETANEG41=PRMTR[19,1];
BETANEG50=PRMTR[20,1];
BETANEG51=PRMTR[21,1];
BETANEG60=PRMTR[22,1];
BETANEG61=PRMTR[23,1];
BETANEG70=PRMTR[24,1];
BETANEG71=PRMTR[25,1];
 
BETAPOS10=PRMTR[26,1];
BETAPOS11=PRMTR[27,1];
BETAPOS20=PRMTR[28,1];
BETAPOS21=PRMTR[29,1];
BETAPOS30=PRMTR[30,1];
BETAPOS31=PRMTR[31,1];
BETAPOS40=PRMTR[32,1];
BETAPOS41=PRMTR[33,1];
BETAPOS50=PRMTR[34,1];
BETAPOS51=PRMTR[35,1];
BETAPOS60=PRMTR[36,1];
BETAPOS61=PRMTR[37,1];
BETAPOS70=PRMTR[38,1];
BETAPOS71=PRMTR[39,1];
 
VAR=PRMTR[40,1]^2;
 
MU_MAT=ST_MAT*MU1 + (ONES(DMNSION,NO_ST)-ST_MAT)*MU0;
 
BETANEG1_MAT=ST_MAT*BETANEG11 + (ONES(DMNSION,NO_ST)-ST_MAT)*BETANEG10;
BETANEG2_MAT=ST_MAT*BETANEG21 + (ONES(DMNSION,NO_ST)-ST_MAT)*BETANEG20;
BETANEG3_MAT=ST_MAT*BETANEG31 + (ONES(DMNSION,NO_ST)-ST_MAT)*BETANEG30;
BETANEG4_MAT=ST_MAT*BETANEG41 + (ONES(DMNSION,NO_ST)-ST_MAT)*BETANEG40;
BETANEG5_MAT=ST_MAT*BETANEG51 + (ONES(DMNSION,NO_ST)-ST_MAT)*BETANEG50;
BETANEG6_MAT=ST_MAT*BETANEG61 + (ONES(DMNSION,NO_ST)-ST_MAT)*BETANEG60;
BETANEG7_MAT=ST_MAT*BETANEG71 + (ONES(DMNSION,NO_ST)-ST_MAT)*BETANEG70;
 
BETAPOS1_MAT=ST_MAT*BETAPOS11 + (ONES(DMNSION,NO_ST)-ST_MAT)*BETAPOS10;
BETAPOS2_MAT=ST_MAT*BETAPOS21 + (ONES(DMNSION,NO_ST)-ST_MAT)*BETAPOS20;
BETAPOS3_MAT=ST_MAT*BETAPOS31 + (ONES(DMNSION,NO_ST)-ST_MAT)*BETAPOS30;
BETAPOS4_MAT=ST_MAT*BETAPOS41 + (ONES(DMNSION,NO_ST)-ST_MAT)*BETAPOS40;
BETAPOS5_MAT=ST_MAT*BETAPOS51 + (ONES(DMNSION,NO_ST)-ST_MAT)*BETAPOS50;
BETAPOS6_MAT=ST_MAT*BETAPOS61 + (ONES(DMNSION,NO_ST)-ST_MAT)*BETAPOS60;
BETAPOS7_MAT=ST_MAT*BETAPOS71 + (ONES(DMNSION,NO_ST)-ST_MAT)*BETAPOS70;
 
/* Matriz de prob. de transição */
 
PR_TR=(QPR~ (1-PPR))|
((1-QPR)~ PPR);
 
/* Inicialização do filtro */
 
A = (eye(2)-pr_tr)|ones(1,2);
EN=(0|0|1);
PROB__T = INV(A'A)*A'EN; //PR[S_t=0]|PR[S_t=1]
 
PR_TRF0=VEC(PR_TR);
PR_TRF1=PR_TRF0|PR_TRF0;
PR_TRF2=PR_TRF1|PR_TRF1;
PR_TRF3=PR_TRF2|PR_TRF2;
PR_TRF4=PR_TRF3|PR_TRF3;
PR_TRF5 =PR_TRF4|PR_TRF4;
PR_TRF =PR_TRF5|PR_TRF5;
 
PROB__T=VECR(PROB__T~PROB__T).*PR_TRF0;
 
PROB__T=VECR(PROB__T~PROB__T).*PR_TRF1;
 
PROB__T=VECR(PROB__T~PROB__T).*PR_TRF2;
 
PROB__T=VECR(PROB__T~PROB__T).*PR_TRF3;
 
PROB__T=VECR(PROB__T~PROB__T).*PR_TRF4;
 
PROB__T=VECR(PROB__T~PROB__T).*PR_TRF5;
 
PROB__=VECR(PROB__T~PROB__T);
 
LIKV=0.0;
J_ITER=1;
DO UNTIL J_ITER>T;
 
F_CAST1=(YY[J_ITER,1]-X_MAT[J_ITER,.]*PHI)*ONES(DMNSION,1)
-(MU_MAT[.,8] - MU_MAT[.,7]*PHI[1,1] - MU_MAT[.,6]*PHI[2,1] - MU_MAT[.,5]*PHI[3,1]
- MU_MAT[.,4]*PHI[4,1] - MU_MAT[.,3]*PHI[5,1]- MU_MAT[.,2]*PHI[6,1] - MU_MAT[.,1]*PHI[7,1])
-Z_MAT[J_ITER,1]*BETANEG1_MAT[.,7] - Z_MAT[J_ITER,2]*BETANEG2_MAT[.,6]
-Z_MAT[J_ITER,3]*BETANEG3_MAT[.,5] - Z_MAT[J_ITER,4]*BETANEG4_MAT[.,4]
-Z_MAT[J_ITER,5]*BETANEG5_MAT[.,3] - Z_MAT[J_ITER,6]*BETANEG6_MAT[.,2]
-Z_MAT[J_ITER,7]*BETANEG7_MAT[.,1]
-Z_MAT[J_ITER,8]*BETAPOS1_MAT[.,7] - Z_MAT[J_ITER,9]*BETAPOS2_MAT[.,6]
-Z_MAT[J_ITER,10]*BETAPOS3_MAT[.,5] - Z_MAT[J_ITER,11]*BETAPOS4_MAT[.,4]
-Z_MAT[J_ITER,12]*BETAPOS5_MAT[.,3] - Z_MAT[J_ITER,13]*BETAPOS6_MAT[.,2]
-Z_MAT[J_ITER,14]*BETAPOS7_MAT[.,1];
 
VAR_L=VAR*ONES(DMNSION,1);
 
PROB_DD=PR_TRF .* PROB__;
 
PR_VL=(1./SQRT(2.*PI.*VAR_L)).*EXP(-0.5*F_CAST1.*F_CAST1./VAR_L).*PROB_DD;
 
PR_VAL=SUMC(PR_VL);
 
LIK=-1*LN(PR_VAL);
 
PRO_=PR_VL/PR_VAL;
 
PROB__T=PRO_[1:DMNSION/2,1]+PRO_[DMNSION/2+1:DMNSION,1];
 
PROB__=VECR(PROB__T~PROB__T);
 
LIKV = LIKV+LIK;
J_ITER = J_ITER+1;
ENDO;
 
LOCATE 2,35;
//"LIKV=";;LIKV;
RETP(LIKV);
ENDP;
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
PROC (3) = FILTER(PRMTR1);
local prmtr, ppr,qpr,prob__0,prob__1,QQ, lik, j_iter, pr__0_l,pr__1_l,
F_cast, var_L,pr_vl, pr_val,likv,phi,PSIC,PSIX,
vecp,st,st1,st2,st3,st4,ST5,ST6,ST7,ST8,ST9,ST10,ST11,ST12,ST13,
pr_tr,pr_trf1,pr_trf,prob__t,prob__,pro_,j,psi1,psi2,var_c,
delta0,DELTA1,MU0,MU1,st_k,st_k1,st_k2,st_k3,st_k4,
f_cast1,f_cast2,PR_VL1,PR_VL2,pr_trf7,pr_trf0,
PR_TRF2,PR_TRF3,PR_TRF4,PR_TRF6,psiL,PR_STL0,PR_STT0,
TMPRY1,TMPRY2,SM_PRL,TMP_P0,SM_PR0,JJJ,MU_MAT,D_MAT,FLT_PRN,
F1,F2,TMP00,TMP0,SM_PR00,prob_dd,VAR,TMP,OUT_MAT,A,EN,BETANEG10,BETANEG11,
BETANEG20,BETANEG21,BETANEG30,BETANEG31,BETANEG40,BETANEG41,BETANEG50,BETANEG51,BETANEG60,BETANEG61,BETANEG70,BETANEG71,BETANEG80,BETANEG81,BETANEG90,BETANEG91,
BETANEG1_MAT,BETANEG2_MAT,BETANEG3_MAT,BETANEG4_MAT,BETANEG5_MAT,BETANEG6_MAT,BETANEG7_MAT,BETANEG8_MAT,BETANEG9_MAT,
BETAPOS10,BETAPOS11,BETAPOS20,BETAPOS21,BETAPOS30,BETAPOS31,BETAPOS40,BETAPOS41,BETAPOS50,BETAPOS51,BETAPOS60,BETAPOS61,
BETAPOS70,BETAPOS71,BETAPOS80,BETAPOS81,BETAPOS90,BETAPOS91,BETAPOS1_MAT,BETAPOS2_MAT,BETAPOS3_MAT,BETAPOS4_MAT,BETAPOS5_MAT,BETAPOS6_MAT,BETAPOS7_MAT,BETAPOS8_MAT,BETAPOS9_MAT, PR_TRF5;
 
PRMTR=TRANS(PRMTR1);
 
LOCATE 15,1; PRMTR';
 
PPR=PRMTR[1,1]; //Pr[St=1/St-1=1] St=1 indica expansao
QPR=PRMTR[2,1]; //Pr[St=0/St-1=0] St=0 indica recessao
 
MU0=PRMTR[3,1]; // tx. cresc. recessao
MU1=PRMTR[4,1]; // tx. cresc. expansao
PHI=PRMTR[5:11,1];
 
BETANEG10=PRMTR[12,1];
BETANEG11=PRMTR[13,1];
BETANEG20=PRMTR[14,1];
BETANEG21=PRMTR[15,1];
BETANEG30=PRMTR[16,1];
BETANEG31=PRMTR[17,1];
BETANEG40=PRMTR[18,1];
BETANEG41=PRMTR[19,1];
BETANEG50=PRMTR[20,1];
BETANEG51=PRMTR[21,1];
BETANEG60=PRMTR[22,1];
BETANEG61=PRMTR[23,1];
BETANEG70=PRMTR[24,1];
BETANEG71=PRMTR[25,1];
 
BETAPOS10=PRMTR[26,1];
BETAPOS11=PRMTR[27,1];
BETAPOS20=PRMTR[28,1];
BETAPOS21=PRMTR[29,1];
BETAPOS30=PRMTR[30,1];
BETAPOS31=PRMTR[31,1];
BETAPOS40=PRMTR[32,1];
BETAPOS41=PRMTR[33,1];
BETAPOS50=PRMTR[34,1];
BETAPOS51=PRMTR[35,1];
BETAPOS60=PRMTR[36,1];
BETAPOS61=PRMTR[37,1];
BETAPOS70=PRMTR[38,1];
BETAPOS71=PRMTR[39,1];
 
VAR=PRMTR[40,1]^2;
 
MU_MAT=ST_MAT*MU1 + (ONES(DMNSION,NO_ST)-ST_MAT)*MU0;
 
BETANEG1_MAT=ST_MAT*BETANEG11 + (ONES(DMNSION,NO_ST)-ST_MAT)*BETANEG10;
BETANEG2_MAT=ST_MAT*BETANEG21 + (ONES(DMNSION,NO_ST)-ST_MAT)*BETANEG20;
BETANEG3_MAT=ST_MAT*BETANEG31 + (ONES(DMNSION,NO_ST)-ST_MAT)*BETANEG30;
BETANEG4_MAT=ST_MAT*BETANEG41 + (ONES(DMNSION,NO_ST)-ST_MAT)*BETANEG40;
BETANEG5_MAT=ST_MAT*BETANEG51 + (ONES(DMNSION,NO_ST)-ST_MAT)*BETANEG50;
BETANEG6_MAT=ST_MAT*BETANEG61 + (ONES(DMNSION,NO_ST)-ST_MAT)*BETANEG60;
BETANEG7_MAT=ST_MAT*BETANEG71 + (ONES(DMNSION,NO_ST)-ST_MAT)*BETANEG70;
 
BETAPOS1_MAT=ST_MAT*BETAPOS11 + (ONES(DMNSION,NO_ST)-ST_MAT)*BETAPOS10;
BETAPOS2_MAT=ST_MAT*BETAPOS21 + (ONES(DMNSION,NO_ST)-ST_MAT)*BETAPOS20;
BETAPOS3_MAT=ST_MAT*BETAPOS31 + (ONES(DMNSION,NO_ST)-ST_MAT)*BETAPOS30;
BETAPOS4_MAT=ST_MAT*BETAPOS41 + (ONES(DMNSION,NO_ST)-ST_MAT)*BETAPOS40;
BETAPOS5_MAT=ST_MAT*BETAPOS51 + (ONES(DMNSION,NO_ST)-ST_MAT)*BETAPOS50;
BETAPOS6_MAT=ST_MAT*BETAPOS61 + (ONES(DMNSION,NO_ST)-ST_MAT)*BETAPOS60;
BETAPOS7_MAT=ST_MAT*BETAPOS71 + (ONES(DMNSION,NO_ST)-ST_MAT)*BETAPOS70;
 
PR_TR=(QPR~ (1-PPR))|
((1-QPR)~ PPR);
 
A = (eye(2)-pr_tr)|ones(1,2);
EN=(0|0|1);
PROB__T = INV(A'A)*A'EN; //PR[S_t=0]|PR[S_t=1]
 
PR_TRF0=VEC(PR_TR);
PR_TRF1=PR_TRF0|PR_TRF0;
PR_TRF2=PR_TRF1|PR_TRF1;
PR_TRF3=PR_TRF2|PR_TRF2;
PR_TRF4=PR_TRF3|PR_TRF3;
PR_TRF5 =PR_TRF4|PR_TRF4;
PR_TRF =PR_TRF5|PR_TRF5;
 
PROB__T=VECR(PROB__T~PROB__T).*PR_TRF0;
 
PROB__T=VECR(PROB__T~PROB__T).*PR_TRF1;
 
PROB__T=VECR(PROB__T~PROB__T).*PR_TRF2;
 
PROB__T=VECR(PROB__T~PROB__T).*PR_TRF3;
 
PROB__T=VECR(PROB__T~PROB__T).*PR_TRF4;
 
PROB__T=VECR(PROB__T~PROB__T).*PR_TRF5;
 
PROB__=VECR(PROB__T~PROB__T);
 
PR_STT0=ZEROS(T,1); // salva Pr[S_t=0|Y_{t}
PR_STL0=ZEROS(T,1); // salva Pr[S_t=0|Y_{t-1}
OUT_MAT=ZEROS(T,1);
 
LIKV=0.0;
J_ITER=1;
DO UNTIL J_ITER>T;
 
F_CAST1=(YY[J_ITER,1]-X_MAT[J_ITER,.]*PHI)*ONES(DMNSION,1)
-(MU_MAT[.,8] - MU_MAT[.,7]*PHI[1,1] - MU_MAT[.,6]*PHI[2,1] - MU_MAT[.,5]*PHI[3,1]
- MU_MAT[.,4]*PHI[4,1] - MU_MAT[.,3]*PHI[5,1]- MU_MAT[.,2]*PHI[6,1] - MU_MAT[.,1]*PHI[7,1])
-Z_MAT[J_ITER,1]*BETANEG1_MAT[.,7] - Z_MAT[J_ITER,2]*BETANEG2_MAT[.,6]
-Z_MAT[J_ITER,3]*BETANEG3_MAT[.,5] - Z_MAT[J_ITER,4]*BETANEG4_MAT[.,4]
-Z_MAT[J_ITER,5]*BETANEG5_MAT[.,3] - Z_MAT[J_ITER,6]*BETANEG6_MAT[.,2]
-Z_MAT[J_ITER,7]*BETANEG7_MAT[.,1]
-Z_MAT[J_ITER,8]*BETAPOS1_MAT[.,7] - Z_MAT[J_ITER,9]*BETAPOS2_MAT[.,6]
-Z_MAT[J_ITER,10]*BETAPOS3_MAT[.,5] - Z_MAT[J_ITER,11]*BETAPOS4_MAT[.,4]
-Z_MAT[J_ITER,12]*BETAPOS5_MAT[.,3] - Z_MAT[J_ITER,13]*BETAPOS6_MAT[.,2]
-Z_MAT[J_ITER,14]*BETAPOS7_MAT[.,1];
 
VAR_L=VAR*ONES(DMNSION,1);
 
PROB_DD=PR_TRF .* PROB__;
 
//-----------------------------------------------------
TMP=PROB_DD;
TMP=TMP[1:128]+TMP[129:256];
TMP=TMP[1:64]+TMP[65:128];
TMP=TMP[1:32]+TMP[33:64];
TMP=TMP[1:16]+TMP[17:32];
TMP=TMP[1:8]+TMP[9:16];
TMP=TMP[1:4]+TMP[5:8];
TMP=TMP[1:2]+TMP[3:4];
 
out_mat[j_iter,.]=(YY[J_ITER,1]-SUMC(PROB_DD.*F_CAST1));
 
PR_STL0[J_ITER,1]=TMP[1,1]; //Pr[S_t=0|Y_t]
//------------------------------------------------------
 
PR_VL=(1./SQRT(2.*PI.*VAR_L)).*EXP(-0.5*F_CAST1.*F_CAST1./VAR_L).*PROB_DD;
 
PR_VAL=SUMC(PR_VL); //f(y_t| Y_{t-1})
LIK=-1*LN(PR_VAL);
 
PRO_=PR_VL/PR_VAL;
 
//-------------------------------------------------------
TMP=PRO_;
TMP=TMP[1:128]+TMP[129:256];
TMP=TMP[1:64]+TMP[65:128];
TMP=TMP[1:32]+TMP[33:64];
TMP=TMP[1:16]+TMP[17:32];
TMP=TMP[1:8]+TMP[9:16];
TMP=TMP[1:4]+TMP[5:8];
TMP=TMP[1:2]+TMP[3:4];
 
PR_STT0[J_ITER,1]=TMP[1,1]; //Pr[S_t=0|Y_t]
//-------------------------------------------------------
 
PROB__T=PRO_[1:DMNSION/2,1]+PRO_[DMNSION/2+1:DMNSION,1];
 
PROB__=VECR(PROB__T~PROB__T);
 
J_ITER = J_ITER+1;
ENDO;
 
RETP(PR_STT0,PR_STL0,OUT_MAT);
ENDP;
//======================================================================
//======================================================================
 
PROC SMOOTH(pr_tt0,pr_tl0);
 
//pr_TT0 contains Pr[S_t|Y_t]
//pr_TL0 contains Pr[S_t|Y_{t-1}]
 
local ppr, qpr, pr_sm0,pr_sm1, j_iter,pr_sm00,pr_sm01,pr_sm10,pr_sm11,
pr_tt1,pr_tl1;
 
PPR=PRM_fnl[1,1]; //Pr[St=1/St-1=1]
QPR=PRM_fnl[2,1]; //Pr[St=0/St-1=0]
 
pr_tt1=1-pr_tt0;
pr_tl1=1-pr_tl0;
 
pr_sm0=pr_tt0;
pr_sm1=pr_tt1;
 
j_iter=T-1;
do until j_iter < 1;
 
pr_sm00=pr_sm0[j_iter+1,1]*qpr* pr_tt0[j_iter,1]/ pr_tl0[j_iter+1,1];
 
pr_sm01=pr_sm1[j_iter+1,1]*(1-qpr)*pr_tt0[j_iter,1]/ pr_tl1[j_iter+1,1];
 
pr_sm10=pr_sm0[j_iter+1,1]*(1-ppr)*pr_tt1[j_iter,1]/ pr_tl0[j_iter+1,1];
 
pr_sm11=pr_sm1[j_iter+1,1]*ppr* pr_tt1[j_iter,1]/ pr_tl1[j_iter+1,1];
 
pr_sm0[j_iter,1]=pr_sm00+pr_sm01;
pr_sm1[j_iter,1]=pr_sm10+pr_sm11;
 
j_iter=j_iter -1;
endo;
 
retp(pr_sm0); //Procedimento gera Pr[S_t=0|Y_T]
endp;
 
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
proc TRANS(c0); // restrição sobre valores das prob. de transição
local c1,m,u,d1,d2,d3,d4,d5,d6;
 
c1=c0;
 
c1[1:2,.]=exp(c0[1:2,.])./ (1+exp(c0[1:2,.]));
 
retp(c1);
endp;