cmlmt calculation error

Just now, I debugged my program and found that the value of the matrixs become a dot after a number of looping. I dont know why and how to correct it. I post my code here, its about 250th loop after which error occurs.

p.s. Could someone tell me how to attach files in this forum?

Thanks for your patience.

---

library cmlmt;

#include cmlmt.sdf;

// load data
r1day1 = xlsReadM("SHIBOR.xlsx","b2:b2059",1,0);
r1day2 = xlsReadM("SHIBOR.xlsx","b3:b2060",1,0);

struct DS d0;
// create a default data structure
d0 = reshape(dsCreate,2,1);
d0[1].datamatrix = r1day1;
d0[2].datamatrix = r1day2;

struct PV p0;
// creates a default parameter structure
p0 = pvCreate;
p0 = pvPackm(p0,0|1,"meanR1",1|1);
// mean parameters in regime 1, arranged as a1,b1
p0 = pvPackm(p0,0|1,"meanR2",1|1);
// mean parameters in regime 2, arranged as a2,b2
p0 = pvPackm(p0,0|.3|.4|1|1,"varR1",1|1|1|1|1);
/* variance parameters in regime1, arranged as
omega1,alpha1,beta1,mu1,v1 */
p0 = pvPackm(p0,0|.4|.3|1|1,"varR2",1|1|1|1|1);
/* variance parameters in regime2, arranged as
omega2,alpha2,beta2,mu2,v2 */
p0 = pvPackm(p0,1|1,"TP",1|1);
// parameters for transition matrix, e1,e2
// p11 = 1/(1+exp(e1)), p22 = 1/(1+exp(e2))
// p0 = pvPackm(p0,1|1,"hstart",1|1);
// parameters for starting value of h(1,0)&h(2,0)
p0 = pvPackm(p0,0|1,"hf1",1|1);
// parameters for calculating f(delata) of regime 1
p0 = pvPackm(p0,0|1,"hf2",1|1);
// parameters for calculating f(delata) of regime 2

struct cmlmtControl c0;
// creates a default structure
c0 = cmlmtControlCreate;
c0.Bounds = { -50 50, -50 50,
-50 50, -50 50,
-10 10, .0001 50, .0001 50, .0001 50, .0001 50,
-10 10, .0001 50, .0001 50, .0001 50, .0001 50,
-50 50, -50 50,
// .0001 10, .0001 10,
-50 50, -50 50,
-50 50, -50 50 };
c0.NumObs = 2058;

// prints the results
call cmlmtprt(cmlmt(&rs_garch,p0,d0,c0));

proc rs_garch(struct pv p, struct ds d, ind);
local meanR1,meanR2,varR1,varR2,TP,hstart,hf1,hf2,tm,
p11,p22,p1,rowN,ex1,ex2,r1dayF,r1dayL,
fvector,pvector,pFvector,pIvector,hvector,hIvector,dIvector,tfvector;
struct modelResults mm;

meanR1 = pvUnpack(p,"meanR1");
/* a1 = meanR1[1];
b1 = meanR1[2]; */
meanR2 = pvUnpack(p,"meanR2");
/* a2 = meanR2[1];
b2 = meanR2[2]; */
varR1 = pvUnpack(p,"varR1");
/* omega1 = varR1[1];
alpha1 = varR1[2];
beta1 = varR1[3];
mu1 = varR1[4];
v1 = varR1[5]; */
varR2 = pvUnpack(p,"varR2");
/* omega2 = varR2[1];
alpha2 = varR2[2];
beta2 = varR2[3];
mu2 = varR2[4];
v2 = varR2[5]; */
TP = pvUnpack(p,"TP");
/* e1 = TP[1];
e2 = TP[2]; */
// hstart = pvUnpack(p,"hstart");
/* h1_0 = hstart[1];
h2_0 = hstart[2]; */
hf1 = pvUnpack(p,"hf1");
/* c1 = hf1[1];
d1 = hf1[2]; */
hf2 = pvUnpack(p,"hf2");
/* c2 = hf2[1];
d2 = hf2[2]; */

p11 = 1/(1+exp(TP[1]));
P22 = 1/(1+exp(TP[2]));
p1 = (1-p22)/(2-p11-p22); // invariant probability of being in regime 1
tm = p11 ~ 1-p11 | 1-p22 ~ p22;// transition matrix

r1dayL = d[1].datamatrix;
r1dayF = d[2].datamatrix;

rowN = rows(r1dayL);

if ind[1];

fvector = zeros(rowN,2); // a matrix consists of [f(r(t+1)|r(t~),s(t)=1),f(r(t+1)|r(t~),s(t)=2)]
pvector = zeros(rowN,2); // a vector consists of [p(s(t)=1|r(t~)),p(s(t)=2|r(t~))]
pIvector = zeros(rowN,2); // a matrix consists of [p1(1,t),p2(1,t)]
pFvector = zeros(rowN,2); // a matrix consists of [p(s(t-1)=1|r(t~)),p(s(t-1)=2|r(t~))]
hvector = zeros(rowN,2); // a matrix consists of [h(1,t),h(2,t)]
hIvector = zeros(rowN,2); // a matrix consists of [hi(1,t),hi(2,t)]
dIvector = zeros(rowN,2); // a matrix consists of [delta_i(1,t),delta_i(2,t)]
tfvector = zeros(rowN,1); // a vector consists of f(r(t+1)|r(t~))

ex1 = r1dayF-(meanR1[1]+meanR1[2]*r1dayL);
ex2 = r1dayF-(meanR2[1]+meanR2[2]*r1dayL);

hvector[1,1] = sumc(ex1^2)/rowN;
hvector[1,2] = sumc(ex2^2)/rowN;

fvector[1,1] = exp(-ex1[1]^2/(2*hvector[1,1]))/(2*pi*hvector[1,1])^.5;
fvector[1,2] = exp(-ex2[1]^2/(2*hvector[1,2]))/(2*pi*hvector[1,2])^.5;

pvector[1,.] = p1~(1-p1);
pFvector[1,.] = pvector[1,.].*fvector[1,.]./sumr(pvector[1,.].*fvector[1,.]);
pIvector[1,1] = p11*pvector[1,1]/(p11*pvector[1,1]+(1-p22)*pvector[1,2]);
pIvector[1,2] = (1-p11)*pvector[1,1]/((1-p11)*pvector[1,1]+p22*pvector[1,2]);
// hvector[1,.] = hstart';
hIvector[1,1] = pIvector[1,1]*hvector[1,1]+(1-pIvector[1,1])*hvector[1,2]
+pIvector[1,1]*(1-pIvector[1,1])*(ex1[1]-ex2[1])^2;
hIvector[1,2] = pIvector[1,2]*hvector[1,1]+(1-pIvector[1,2])*hvector[1,2]
+pIvector[1,2]*(1-pIvector[1,2])*(ex1[1]-ex2[1])^2;
dIvector[1,1] = pIvector[1,1]*ex1[1]/hvector[1,1]^.5+(1-pIvector[1,1])*ex2[1]/hvector[1,2]^.5;
dIvector[1,2] = pIvector[1,2]*ex1[1]/hvector[1,1]^.5+(1-pIvector[1,2])*ex2[1]/hvector[1,2]^.5;

tfvector[1] = sumr(fvector[1,.].*pvector[1,.]);

for i (2,rowN,1);
hvector[i,1] = (varR1[1]+varR1[2]*(hIvector[i-1,1])^(varR1[4]/2) * (abs(dIvector[i-1,1]-hf1[2])-hf1[1]*(dIvector[i-1,1]-hf1[2]))^varR1[5]
+varR1[3]*(hIvector[i-1,1])^(varR1[4]/2))^(2/varR1[4]);
hvector[i,2] = (varR2[1]+varR2[2]*(hIvector[i-1,2])^(varR2[4]/2) * (abs(dIvector[i-1,2]-hf2[2])-hf2[1]*(dIvector[i-1,2]-hf2[2]))^varR2[5]
+varR2[3]*(hIvector[i-1,2])^(varR2[4]/2))^(2/varR2[4]);

fvector[i,1] = exp(-ex1[i]^2/(2*hvector[i,1]))/(2*pi*hvector[i,1])^.5;
fvector[i,2] = exp(-ex2[i]^2/(2*hvector[i,2]))/(2*pi*hvector[i,2])^.5;

pvector[i,.] = pFvector[i-1,.]*tm;
// pvector[i,1] = p11*pFvector[i-1,1]+(1-p22)*pFvector[i-1,2];
// pvector[i,2] = p22*pFvector[i-1,2]+(1-p11)*pFvector[i-1,1];

pFvector[i,.] = pvector[i,.].*fvector[i,.]./sumr(pvector[i,.].*fvector[i,.]);
// pFvector[i,1] = pvector[i,1]*fvector[i,1]/(pvector[i,1]*fvector[i,1]+pvector[i,2]*fvector[i,2]);
// pFvector[i,2] = pvector[i,2]*fvector[i,2]/(pvector[i,1]*fvector[i,1]+pvector[i,2]*fvector[i,2]);

pIvector[i,1] = p11*pvector[i,1]/(p11*pvector[i,1]+(1-p22)*pvector[i,2]);
pIvector[i,2] = (1-p11)*pvector[i,1]/((1-p11)*pvector[i,1]+p22*pvector[i,2]);

hIvector[i,1] = pIvector[i,1]*hvector[i,1]+(1-pIvector[i,1])*hvector[i,2]
+pIvector[i,1]*(1-pIvector[i,1])*(ex1[1]-ex2[1])^2;
hIvector[i,2] = pIvector[i,2]*hvector[i,1]+(1-pIvector[i,2])*hvector[i,2]
+pIvector[i,2]*(1-pIvector[i,2])*(ex1[1]-ex2[1])^2;
dIvector[i,1] = pIvector[i,1]*ex1[1]/hvector[i,1]^.5+(1-pIvector[i,1])*ex2[1]/hvector[i,2]^.5;
dIvector[i,2] = pIvector[i,2]*ex1[1]/hvector[i,1]^.5+(1-pIvector[i,2])*ex2[1]/hvector[i,2]^.5;

tfvector[i] = sumr(fvector[i,.].*pvector[i,.]);

endfor;

mm.function = sumc(ln(tfvector));
endif;
retp(mm);
endp;