/* construction of fractional lagged terms */
new; cls;
y=rndn(100,1);
d=0.4;
p=3; /* lag order */
print vlag(y,d,p);



proc vlag(y,d,p);
   local n,k,lagy,i,dy,y_d,j;
   n=rows(y);
   k=cols(y);
   lagy=miss(zeros(n,k*p),0);
   i=1;
   do while i<=k;
      dy=fracdiff(y[.,i],d); 
      dy=dy[2:n];
      y_d=y[2:n,i]-dy;
      lagy[2:n,i]=y_d;
      i=i+1;
   endo;
   j=2;
   do while j<=p;
      i=1;
      do while i<=k;
         dy=fracdiff(lagy[j:n,k*(j-2)+i],d); 
         dy=dy[2:n-j+1];
         lagy[j+1:n,k*(j-1)+i]=lagy[j+1:n,k*(j-2)+i]-dy;
         i=i+1;
      endo;
      j=j+1;
   endo;
   lagy=lagy[p+1:n,.];
   retp(lagy);
endp;

proc fracdiff(y,d);
   local n,t,y1,x;
   n=rows(y);
   t=seqa(1,1,n-1);
   y1=1|cumprodc((t-1-d)./t);
   x=conv(y,y1,0,n);
   retp(x);
endp;