/* ****************************************************************************
** GAUSS code
** File: tarfit.src                                                       
** Coded by: R. Schlittgen
**
** Reference:
** R. Schlittgen (1997)
**    Fitting of threshold models for time series
**    Discussion paper 97-13, 
**    Department of Economics, University of California, San Diego
**    Available at: http://papers.ssrn.com/sol3/papers.cfm?abstract_id=52401
**
** Function for modeling threshold autoregressive time series
** Format                                 Purpose                    Line
** ------------------------------------------------------------------------
** {par1,par2,r,d,res} = tarfit(x,p,d);     Fits TAR(p1,p2) model       32
*/
/* TARFIT
**
** Purpose: Fits SETAR[p1,p2,2] model.
**
** Reference:
**
**    Format: {par1,par2,c,d,aic,sig} = tarfit(x,p,d);
**
**    Input:  x    : (n,1)-vector, the time series
**            p    : maximum of orders of the autoregressive regimes
**            d    : maximum of delay
**    Output: par1 : (p1+1,1)-vector, parameters of first AR-regime
**            par2 : (p2+1,1)-vector, parameters of second AR-regime
**            c    : scalar, threshold parameter
**            d    : scalar, delay parameter
**            aic  : minimum of modified AIC
**            sig  : pooled estimate of error variance
******************************************************************** */
proc (6) = tarfit(x,p,d);
   local a,aic,aic1,aic2,c,dd,dmat,ind,m,mat,mat1,mini;
   local n,ns,par1,par2,p1,p2,r,res,R1,R2,s,sig,sig1,sig2,s1,s2;

   dd = 0;
   dmat = zeros(d,4);

   n = rows(x);
   if n==1;
      x=x';
      n=rows(x);
   endif;

   mat1 = trimr(rev(reshape(x,p+1,n+1))',0,p+1);
   mat1 = ones(n-p,1)~mat1[.,2:p+1]~mat1[.,1];
   s = seqa(1,1,p+1);

   do while dd<d;

      dd = dd+1;

      if dd <= p;
         mat =  sortc(mat1~x[p-dd+1:n-dd],p+3);
      else;
         mat = sortc(trimr(mat1,dd-p,0)~x[1:n-dd],p+3);
      endif;

      m    = rows(mat);
      c    = mat[.,p+3];
      mat  = mat[.,1:p+2];
      mini = p+4;
      r    = mini;
      ns   = r+1+p-s;

      R1   = qr(mat[1:r,.]);
      sig1 = R1[p+2,p+2]^2;
      a    = s~ns~(ln(rev(cumsumc(R1[p+2:2,p+2]^2))./ns) + 2*s./ns);
      a    = a[minindc(a[.,3]).,];
      aic1 = c[r]~a;

      R2   = qr(mat[m+1-r:m,.]);
      sig2 = R2[p+2,p+2]^2;
      a    = s~ns~(ln(rev(cumsumc(R2[p+2:2,p+2]^2))./ns) + 2*s./ns);
      a    = a[minindc(a[.,3]).,];
      aic2 = a;

      do while r < m-mini;
        r  = r+1;
        ns = ns+1;

        R1          = qr(R1[1:p+1,.]|mat[r,.]);
        sig1        = sig1 + R1[p+2,p+2]^2;
        R1[p+2,p+2] = sqrt(sig1);
        a           = s~ns~(ln(rev(cumsumc(R1[p+2:2,p+2]^2))./ns) + 2*s./ns);
        a           = a[minindc(a[.,3]).,];
        aic1        = aic1|(c[r]~a);

        R2          = qr(R2[1:p+1,.]|mat[m+1-r,.]);
        sig2        = sig2 + R2[p+2,p+2]^2;
        R2[p+2,p+2] = sqrt(sig2);
        a           = s~ns~(ln(rev(cumsumc(R2[p+2:2,p+2]^2))./ns) + 2*s./ns);
        a           = a[minindc(a[.,3]).,];
        aic2        = a|aic2;
      endo;

      aic = aic1[.,1 2]~aic2[.,1];
      aic = aic~((aic1[.,3].*aic1[.,4] +
                  aic2[.,2].*aic2[.,3])./(aic1[.,3]+aic2[.,2]));
      dmat[dd,.] = aic[minindc(aic[.,4]),.];
    endo;

    dd = minindc(dmat[.,4]);
    p1 = dmat[dd,2]-1;
    p2 = dmat[dd,3]-1;
    c  = dmat[dd,1];

    if p1==0;                                          
       mat = ones(n,1)~x;                              
    else;                                              
       mat = trimr(rev(reshape(x,p1+1,n+1))',0,p1+1);
       mat = ones(n-p1,1)~mat[.,2:p1+1]~mat[.,1];
    endif;                                             

    if dd<=p1;
       mat = selif(mat,x[p1-dd+1:n-dd].<=c);
    else;
       mat = selif(trimr(mat,dd-p1,0),x[1:n-dd].<=c);
    endif;

    m    = rows(mat);
    R1   = qr(mat);
    sig1 = R1[p1+2,p1+2]^2;
    aic1 = m*ln(sig1/m) + 2*(p1+1);
    par1 = qrsol(R1[1:p1+1,p1+2],R1[1:p1+1,1:p1+1]);

    if p2==0;                                          
       mat = ones(n,1)~x;                              
    else;                                              
       mat = trimr(rev(reshape(x,p2+1,n+1))',0,p2+1);
       mat = ones(n-p2,1)~mat[.,2:p2+1]~mat[.,1];
    endif;                                             

    if dd<=p2;
       mat = selif(mat,x[p2-dd+1:n-dd].>c);
    else;
       mat = selif(trimr(mat,dd-p2,0),x[1:n-dd].>c);
    endif;

    m    = rows(mat);
    R2   = qr(mat);
    sig2 = R2[p2+2,p2+2]^2;
    aic2 = m*ln(sig2/m) + 2*(p2+1);
    par2 = qrsol(R2[1:p2+1,p2+2],R2[1:p2+1,1:p2+1]);

    aic  = (aic1+aic2)/(n-dd);
    @ sig = (sig1+sig2)/(n-dd-2-p1-p2);@
    sig  = (sig1+sig2)/n;

  retp(par1,par2,c,dd,aic,sig);
endp;
