clear dat1,xx,wei;   @Globals for Pantob@
clear _amftol,_amalp,_ambet,_amgam,_amitmax,__output;   @Globals for Amoeba@

_amftol=0.000001;  _amalp=1;   _ambet=0.5;    @Default values@
_amgam=2.0; _amitmax=200; _ammaxt=1.0d10;

clear __miss;

#include optmum.prc;

proc bopt1(y,xb,nobs,xx);
local f, i,is,it,xbdif;

f=0;
for i (1,nobs,1);
for is (xx[i,1],xx[i,2]-1,1);
for it (is+1,xx[i,2],1);
xbdif=xb[is]-xb[it];
if xbdif.lt -1; 
    xbdif=-1; 
elseif xbdif.gt 1;
    xbdif=1;
endif;
f=f+wei[i]*fct_r(y[is],y[it],xbdif);
endfor;
endfor;
endfor;
retp(-f);
endp;

proc bodpt1(y,xb,x,n,nobs,xx,k);
local df,i,is,it,xbdif;

df=zeros(k,1);
for i (1,nobs,1);
for is (xx[i,1],xx[i,2]-1,1);
for it (is+1,xx[i,2],1);
xbdif=xb[is]-xb[it];
if xbdif.lt -1; 
    xbdif=-1; 
elseif xbdif.gt 1;
    xbdif=1;
endif;
df=df+wei[i]*fct_u(y[is],y[it],xbdif)*(x[is,.]-x[it,.])';
endfor;
endfor;
endfor;
retp(-df);
endp;

proc bodpt1a(y,xb,x,n,nobs,xx,k);
local df,i,is,it,xbdif;

df=zeros(k,nobs);
for i (1,nobs,1);
for is (xx[i,1],xx[i,2]-1,1);
for it (is+1,xx[i,2],1);
xbdif=xb[is]-xb[it];
if xbdif.lt -1; 
    xbdif=-1; 
elseif xbdif.gt 1;
    xbdif=1;
endif;
df[.,i]=df[.,i]+wei[i]*fct_u(y[is],y[it],xbdif)*(x[is,.]-x[it,.])';
endfor;
endfor;
endfor;
retp(-df);
endp;

proc fct_r(y1,y2,d);
local i1,i2,i3,i4,f1,f2,f3,f4,u1,u2;



f1=d+d.*d/2+(y1-1)^2/2;
f2=d*y1;
f3=y1*d-d.*d/2;
f4=y1*y1/2;

i1=d.lt (y1-1);
i2=(d.lt 0).*(d.ge (y1-1));
i3=(d.ge 0).*(d.lt y1);
i4=(d.ge y1);
u1=f1.*i1+f2.*i2+f3.*i3+f4.*i4;


f1=-y2.*y2/2;
f2=d.*d/2+d*y2;
f3=d*y2;
f4=d-d.*d/2-(1-y2)+(1-y2)*(1-y2)/2+(1-y2)*y2;
f4=d-d.*d/2-(1-y2)*(1-y2)/2;

i1=d.lt -y2;
i2=(d.lt 0).*(d.ge -y2);
i3=(d.ge 0).*(d.lt (1-y2));
i4=(d.ge (1-y2));
u2=f1.*i1+f2.*i2+f3.*i3+f4.*i4;


retp(u1-u2);
endp;

proc fct_u(y1,y2,d);
local a1,a2,u1,u2;

a1=(y1-d);
a1=a1.*(a1.gt 0);
a2=y1.lt (1+d);
a2=y1.*a2 + (1+d).*(1-a2);

u1=a1.*(d.gt 0) +a2.*(d.le 0);

a1=(y2+d);
a1=a1.*(a1.gt 0);

a2=y2.lt (1-d);
a2=y2.*a2 + (1-d).*(1-a2);

u2=a2.*(d.gt 0) +a1.*(d.le 0);

retp(u1-u2);
endp;

proc bofu1(b);
/* function evaluation */
local xb,y,nobs,f;
    nobs=rows(xx);
    xb=dat1[.,3:cols(dat1)]*b;
    f=0;
    y=dat1[.,1];
    f=bopt1(y,xb,nobs,xx);

retp(f);
ENDP;

PROC bodfu1(b);
/* first derivatives */
local xb,y,nobs,df,dfi,x,k,n;
    nobs=rows(xx);
    n=rows(dat1);
    x=dat1[.,3:cols(dat1)];
    xb=x*b;
    y=dat1[.,1];
    k=cols(x);
    df=zeros(k,1);  
    df=bodpt1(y,xb,x,n,nobs,xx,k);
retp(df);
endp;

proc bodfu1a(b);
/* first derivatives */
local xb,y,nobs,df,dfi,x,k,n;
    nobs=rows(xx);
    n=rows(dat1);
    x=dat1[.,3:cols(dat1)];
    xb=x*b;
    y=dat1[.,1];
    k=cols(x);
    df=bodpt1a(y,xb,x,n,nobs,xx,k);
retp(df);
endp;



PROC boddfu1(b);
/* second derivatives */
local xb,y,nobs,ga,x,k,n,i1,i2,i3,i4,is,it,xbdif,dx;
    nobs=rows(xx);
    n=rows(dat1);
    x=dat1[.,3:cols(dat1)];
    xb=x*b;
    y=dat1[.,1];
    k=cols(x);
    ga=zeros(k,k);

    for i (1,nobs,1);
    for is (xx[i,1],xx[i,2]-1,1);
    for it (is+1,xx[i,2],1);
    xbdif=xb[is]-xb[it];
    if xbdif.lt -1; 
        xbdif=-1; 
    elseif xbdif.gt 1;
        xbdif=1;
    endif;    
    i1=(-1.lt xbdif).*(xbdif.lt (y[is]-1));
    i2=(0.lt xbdif).*(xbdif.lt y[is]);
    i3=(-y[it].lt xbdif).*(xbdif.lt 0);
    i4=((1-y[it]).lt xbdif).*(xbdif.lt 1);
    dx=(x[is,.]-x[it,.])';
    ga=ga+wei[i]*(i1-i2-i3+i4)*(dx*dx');
    endfor;
    endfor;
    endfor;

retp(-ga);
ENDP;

proc(4)=sortda1(z,i);
local x1,x2,zz,n,ii,wei;
    z=sortc(z,i);
    n=rows(z);
    x1=zeros(n,1);
    x2=zeros(n,1);
    zz=z[.,i];

    if (n.le 0);
    print "ERROR - No individual has more than one observation. The data is not
          a true panel.";
    else;
    ii=1;
    x1[ii]=1; x2[ii]=0;
    for i (2,n,1);
        if zz[i].ne zz[i-1];
            if (x1[ii].eq (i-1)); 
                x1[ii]=i;
            else;
                x2[ii]=i-1;
                ii=ii+1;
                x1[ii]=i;
            endif;
        endif;
    endfor;
    endif;
    x2[ii]=n;
    n=ii;
    wei=1+x2[1:n]-x1[1:n];
    x1=x1[1:n]~x2[1:n];


retp(z,x1,wei,n);
endp;

proc(3)=twoside(x,_ptmet);
      local fh1,varnames,xnames,yname,iname,z,zout,k,n,startp,starty,point,
        ttt,pointval,pstar,ystar,it,ic,ind,ilo,bhat,f,g,g_1,v,vcv,b,n_1,imet,
        ff,b0;

    b0=x[.,1]/x[.,3:cols(x)]; 

    n_1=rows(x);
    x=packr(x);
    x=sortc(x,2);

    format /rd 4,0;


   
    if rows(x)<=1;   print "ERROR - All data has been deleted. ";
                    RETP(0,0);
    else;

    {dat1,xx,wei,n}=sortda1(x,2);
     format /rd 4,0;
    " ";
    " The dataset contains" n "individuals.";
 format /ld 4,0;
   n_1-rows(x) "observations of " n_1 "deleted due to missing data."; 
  


    if _ptmet==0;   @   Use Amoeba. @

        k=cols(dat1)-2;
        startp=(zeros(1,k)|diagrv(eye(k),0.1*abs(b0)+1))+b0' ;
        {pstar,ystar,it,ttt}=amoeba(startp,_amftol,_ammaxt,_amitmax,&bofu1);
        bhat=pstar[1,.]';
        f=ystar[1];

    elseif _ptmet==-1;       @  Use Powell's method.    @

        {bhat,ff,ic,ttt} = powell(b0,_amftol,_ammaxt,_amitmax,&bofu1,&linmin);
        f=ff;

    elseif _ptmet==-2;       @  Use frpmin.             @

        {bhat,ff,ic,ttt} = frprmn(b0,_amftol,_ammaxt,_amitmax,&bofu1,&bodfu1,&linmin);
        f=ff;

    elseif _ptmet==-3;       @  Use dfpmin.             @

        {bhat,ff,ic,ttt} = dfpmin(b0,_amftol,_ammaxt,_amitmax,&bofu1,&bodfu1,&linmin);
        f=ff;

    else ;
        "minimizer not found";
    endif;

    g=bodfu1a(bhat);
    v=g*g';
    g_1=boddfu1(bhat);
    if det(g_1)==0;  "d_1 singular"; endif;
    g_1=inv(g_1);
    vcv=g_1'*v*g_1;
print /flush " ";
retp(bhat,vcv,bofu1(bhat));
endif;
ENDP;