clear dat1,xx,_pttheta,_ptloss,_ptbw;   @Globals for Pantob@
_pttheta=3; _ptloss=1;  _ptbw=0;    _ptmet=1;   @Default values@

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;


PROC bofu1(b);
/* function evaluation */
local xb,y,nobs,f;
    nobs=rows(xx);
    xb=dat1[.,3:cols(dat1)]*b;
    f=0;
    y=dat1[.,1];
    dllcall ibopt1(y,xb,nobs,xx,f,_ptloss,_pttheta);
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);  dfi=df;
    dllcall ibodpt1(y,xb,x,n,nobs,xx,df,dfi,k,_ptloss,_pttheta);
RETP(df');
ENDP;

PROC bodfu1p(b);
/* first derivatives as column */
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);  dfi=df;
    dllcall ibodpt1(y,xb,x,n,nobs,xx,df,dfi,k,_ptloss,_pttheta);
RETP(df);
ENDP;

PROC boddfu1(b);
/* second derivatives */
local xb,y,nobs,ga,x,k,n,ga1;
    nobs=rows(xx);
    n=rows(dat1);
    x=dat1[.,3:cols(dat1)];
    xb=x*b;
    y=dat1[.,1];
    k=cols(x);
    ga=zeros(k,k);ga1=ga;
    dllcall iboddpt1(y,xb,x,n,nobs,xx,ga,ga1,k,_ptloss,_pttheta);
RETP(ga);
ENDP;

PROC bonddfu1(b,h);
/* numrical second derivatives */
local b1,b2,k,i,f1,f2,ga;
k=rows(b);
ga=zeros(k,k);
i=1;
if rows(h)==1;
    do while i<=k;
        b1=b; b1[i]=b1[i]-h;
        b2=b; b2[i]=b2[i]+h;
        ga[i,.]=(bodfu1(b2)-bodfu1(b1))/(2*h);
        i=i+1;
    endo;

elseif rows(h)==k;
    do while i<=k;
        b1=b; b1[i]=b1[i]-h[i];
        b2=b; b2[i]=b2[i]+h[i];
        ga[i,.]=(bodfu1(b2)-bodfu1(b1))/(2*h[i]);
        i=i+1;
    endo;

else;
    "Bandwidth vector and parameter vector do not match.";

endif;

RETP(ga);
ENDP;

PROC bocv1(b,d);
/* Vhat matrix */
local xb,y,nobs,df,x,k,n,v;
    nobs=rows(xx);
    n=rows(dat1);
    x=dat1[.,3:cols(dat1)];
    xb=x*b;
    y=dat1[.,1];
    k=cols(x);
    df=zeros(k,1);
    v=zeros(k,k);
    dllcall ibocalv1(y,xb,x,n,nobs,xx,df,k,_ptloss,_pttheta,v);
RETP(v);
ENDP;

PROC(3)=sortda(z,i);
local x1,x2,zz,n;
    z=sortc(z,i);
    n=rows(z);
    x1=zeros(n,1);
    x2=zeros(n,1);
    zz=z[.,i];
    dllcall ibosort(zz,x1,x2,n);
    if (n.le 0);
    print "ERROR - No individual has more than one observation. The data is not
          a true panel.";
    endif;
    x1=x1[1:n]~x2[1:n];
RETP(z,x1,n);
ENDP;


PROC(9)=pantob1(infile,yvar,i,xvars,b0);
      local x,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;

      open fh1 = ^infile;
      x=readr(fh1,100);
      do until eof(fh1);
      x=x|readr(fh1,100);
      endo;

    varnames=getname(infile);

    fh1=close(fh1);

    xnames=varnames[xvars];
    yname=varnames[yvar];
    iname=varnames[i];

    x=x[.,yvar]~x[.,i]~x[.,xvars];
    n_1=rows(x);
    x=packr(x);
    x=sortc(x,2);


    format /rd 4,0;
    cls;
    print "Data from file " infile;
    print n_1-rows(x) "observations of " n_1
            "deleted due to missing data.";
    format /mb1 /ros 16,8;
    if rows(x)<=1;   print "ERROR - All data has been deleted. ";
                    RETP(0,0);
    else;
    "Sorting Data";
    {dat1,xx,n}=sortda(x,2);

 

    if _ptmet==0;   @   Use Amoeba. @

        k=cols(dat1)-2;

        startp=(zeros(1,k)|diagrv(eye(k),0.5*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,&bodfu1p,&linmin);
    f=ff;

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

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


    else ;
"UUPPPSSS";

    endif;


    b0=bhat;

    if _ptbw==0 and _ptloss==2;
    "Error: Bandwidth must be positive for absolute value loss function";
    "Setting _ptbw=0.125 and continuing."; _ptbw=0.125;
    endif;

    v=bocv1(bhat,0);
    if _ptbw==0; g_1=boddfu1(bhat);
    elseif _ptbw>0; g_1=bonddfu1(bhat,_ptbw);
    endif;

    if det(g_1)==0;  output file=g_1error on;  g_1;  output off; endif;

    g_1=inv(g_1);
    vcv=g_1'*v*g_1;

RETP(bhat,vcv,f,infile,_ptmet,yname,iname,xnames,n);
endif;
ENDP;