new;
cls;

PROC cdfe(x);
local xla;
xla=exp(x*pi/sqrt(3));
RETP(xla./(1+xla));
ENDP;

PROC bosec;
local tm;
tm=time;
RETP(tm[1]*3600+tm[2]*60+tm[3]+tm[4]/100);
ENDP;

PROC(2)=probs(g,b,al);
local p0,p1,d1,d2,d3,c1g1,c1g0,pr2,pr3
         ,ct11,ct10,ct01,ct00,t,it,c2,c3,nseq,jj,c4,pr4,d4;

nseq=8;
t=3;

p0=zeros(nseq,rows(al));
p1=p0;

ct11=zeros(rows(al),t);
ct10=ct11;
ct01=ct11;
ct00=ct11;

it=1;
do while it.le t;
ct11[.,it]=cdfn(g+b[it]+al);
ct01[.,it]=cdfn(b[it]+al);
it=it+1;
endo;
ct10=1-ct11;
ct00=1-ct01;

jj=1;
for d1 (0,1,1);

/*  prob of observed in first period condition on 0th period */

c1g1=(ct11[.,1]*d1+ct10[.,1]*(1-d1));
c1g0=(ct01[.,1]*d1+ct00[.,1]*(1-d1));

/*  prob next period condition on observed this period */

pr2= ct11[.,2]*d1+ct01[.,2]*(1-d1);

for d2 (0,1,1);

c2=pr2*d2+(1-pr2)*(1-d2);
pr3= ct11[.,3]*d2+ct01[.,3]*(1-d2);

for d3 (0,1,1);

c3=pr3*d3+(1-pr3)*(1-d3);

p1[jj,.]=(c1g1 .*c2 .* c3)' ;
p0[jj,.]=(c1g0 .*c2 .* c3)' ;

jj=jj+1;


endfor;
endfor;
endfor;


RETP(p1,p0);
ENDP;


dlibrary c:\research\forgau95\sim_new\release\sim_new.dll;

gamm=seqa(-1,0.01,201);
maxg=gamm*0;
ming=maxg;

for ig (1,rows(gamm),1);
ga=gamm[ig];  b_tru=0.0;
beta=0|b_tru|(2*b_tru);
n_al=101;

al=seqa(-10,0.2,n_al); 
p_al=zeros(n_al,1);
p_al[1]=cdfn((al[1]+al[2])/2);
i_al=2;
do while i_al.lt n_al;
p_al[i_al]=cdfn((al[i_al]+al[i_al+1])/2)-cdfn((al[i_al-1]+al[i_al])/2);
i_al=i_al+1;
endo;
p_al[n_al]=1-sumc(p_al);

al~p_al;
sumc(al'*p_al);; sumc((al.^2)'*p_al)-sumc(al'*p_al)^2;

p_0=.5*ones(n_al,1); 


{p1,p0}=probs(ga,beta,al);
pdat=p1*(p_0.*p_al)+p0*((1-p_0).*p_al);


botim=bosec;

npoi=501;
dpoi=2/(npoi-1);
sss1=seqa(-1,dpoi,npoi);
dga=1;


res=1000+zeros(1,2);
res[1,1]=-1;
"start";
print /flush " ";

kk=1;

for igam (1,rows(sss1),1);
bt=bosec;

g=ga+sss1[igam]*dga;

bbb=0;
be=0|bbb|(2*bbb); 
{p1,p0}=probs(g,be,al);
aa=p1~p0;
a=aa|ones(1,cols(aa));

pdata=pdat|1;
/* these are the constrains */

xm=rows(a); xn=cols(a);
aa=(0~ones(1,cols(a)) )| (pdata~-a)|zeros(1,(cols(a)+1));
aa=aa';
xmp=cols(aa); xnp=rows(aa); xm1=0; xm2=0; xm3=rows(pdata);
xica=0; xizro=zeros(1,xn); xiposv=zeros(xm,1);

eps=1.0d-6;
dllcall isimpl(aa,xm,xn,xmp,xnp,xm1,xm2,xm3,xica,eps); 

if (xica.ne -1); res=res|(xica~g); endif;

kk=kk+1;




endfor;



bt=bosec;

"HEREeeewf timqwe";  

bt-botim;

rr=selif(res[.,2],res[.,1].eq 0);
ming[ig]=minc(rr);
maxg[ig]=maxc(rr);

endfor;

library pgraph;
graphset;
_pltype=6;
_plwidth=2;
_pcolor=1;
_plctrl=0;
xy(gamm,ming~maxg);