ADF test conversion from Matlab to R

Question

Online, I found the matlab code for the ADF test together with the tests for the SADF and DF test. Now, in order to convert this matlab code to R, line 45 still tend to give me an error.

    # if mflag=1 use ADF with constant and without trend
    # if mlag=2 use ADF with constant and trend
    # if mlag=3 use ADF without constant and trend

    #ADF  

    #Set data sample equal to y to start
    y=Data$OIL.P

    #Start ADF test 
    ADF_FL<- function(y,adflag,mflag)
    {
      t1=dim(y,1)-1
      const=matrix(1,t1,1)
      trend=t(seq(1,1,by=t1)) #or trend=t(c(seq(1,1,by=t)1)))

      y1=y(dim(y,1)-t1:dim(y,1)-1)
      dy=y(2:dim(y, 1))-y(1:dim(y,1)-1)
      dy0=dy(dim(dy,1)-t1+1:dim(dy,1))
      x=y1

      if(mflag==1){
      x=c(x,const)
      } else{
      x=c(x, const, trend)
      }

      x1=x

      t2=t1-adflag
      x2=x1(dim(x1,1)-t2+1:dim(x1,1),) #from k+1 to the end (including y1 and x)
      dy01=dy(dim(dy0,1)-t2+1:dim(dy0,1)) #from k+1 to the end (including dy0)

      if(adflag>0){
        j=1
        while(j<=adflag){
          x2=[x2 dy(dim(dy,1)-t2+1-j:dim(dy,1)-j)] #problem in this line
          j=j+1
        }
      }

      beta=((t(x2))*x2)^(-1)*(t(x2)*dy01)
      eps=dy01-x2*beta

      if(mflag==1){
        sign=sqrt(diag(t(eps)*eps/(t2-adflag-2)*(t(x2)*x2)^(-1)))
      } else{
        sign=sqrt(diag(t(eps)*eps/(t2-adflag-3)*(t(x2)*x2)^(-1)))
      }

      tvalue=beta/sign
      show(tvalue)
     }

Original matlab code:

    function [estm]=ADF_FL(y,adflag,mflag)

      t1=size(y,1)-1;
      const=ones(t1,1);
      trend=1:1:t1;trend=trend';

      y1  = y(size(y,1)-t1:size(y,1)-1);
      dy  = y(2:size(y,1)) - y(1:size(y,1)-1);
      dy0 = dy(size(dy,1)-t1+1:size(dy,1));
      x=y1;

      if mflag==1;  x=[x const]; end;
      if mflag==2;  x=[x const trend]; end;
      % if mflag==3;  x=x; end;

      x1=x; 

      t2=t1-adflag;
      x2=x1(size(x1,1)-t2+1:size(x1,1),:);         % from k+1 to the end  (including y1 and x)-@
      dy01=dy0(size(dy0,1)-t2+1:size(dy0,1));      % from k+1 to the end   (including dy0)-@

      if adflag>0;
        j=1; 
        while j<=adflag; 
        x2=[x2 dy(size(dy,1)-t2+1-j:size(dy,1)-j)];     %including k lag variables of dy in x2-@
        j=j+1; 
      end;
    end;

       beta =(x2'*x2)^(-1)*(x2'*dy01);                       % model A-@
       eps  = dy01 - x2*beta;

       if mflag==1; sig = sqrt(diag(eps'*eps/(t2-adflag-2)*(x2'*x2)^(-1))); end;
       if mflag==2; sig = sqrt(diag(eps'*eps/(t2-adflag-3)*(x2'*x2)^(-1))); end;
       if mflag==3; sig = sqrt(diag(eps'*eps/(t2-adflag-1)*(x2'*x2)^(-1))); end;

       tvalue=beta./sig;

       estm=tvalue(1);

Could someone help me in this?

Thanks,

David

Answer 1

I have that full code, coded in R here it is. With this kind of coding it assumes that y is a vector, and that's why the first step is converting it into a matrix. It works perfectly to find the ADF statistic, and i compare it with the fUnitRoot adfTest and it gives the same results.

ADF_FL = function(y,adflag,mflag){
  y = as.matrix(y, nrow=1)
  t1 = nrow(y)-1
  constant = matrix(1,nrow=t1,ncol=1)
  trend = as.matrix(seq(1,t1,by=1),ncol=1)

  y1 = as.matrix(y[(nrow(y)-t1):(nrow(y)-1),],ncol=1)
  dy =as.matrix(y[-1,]- y[1:(nrow(y)-1),],ncol=1)
  dy0 = as.matrix(dy[(nrow(dy)-t1+1):nrow(dy)],ncol=1)
  x = y1

  if(mflag == 1){
    x = cbind(constant,x)
  }
  if(mflag == 2){
    x = cbin(constant,trend,x)
  } else{x = x}

  x1 = x

  t2 = t1 - adflag
  x2 = x1[(nrow(x1)-t2+1):nrow(x1),]       
  dy01 = as.matrix(dy0[(nrow(dy0)-t2+1):nrow(dy0),],ncol=1)

  if(adflag>0){
    j = 1
    while(j<=adflag){
      x2 = cbind(x2,dy[(nrow(dy)-t2+1-j):nrow(dy)-j])
      j = j+1
    }
  }

  beta = solve(t(x2)%*%x2)%*%(t(x2)%*%dy01) ## beta de la regresion de la prueba ADF
  eps = dy01-x2%*%beta ## errores de la regresion 
  if(mflag == 1){sig = sqrt(diag(as.numeric(t(eps)%*%(eps/(t2-adflag-2)))*solve(t(x2)%*%x2)))}
  if(mflag == 2){sig = sqrt(diag(as.numeric(t(eps)%*%(eps/(t2-adflag-3)))*solve(t(x2)%*%x2)))}
  if(mflag == 3){sig = sqrt(diag(as.numeric(t(eps)%*%(eps/(t2-adflag-1)))*solve(t(x2)%*%x2)))}

  tvalue = beta/sig
  return(tvalue[2])
}