HOW TO RUN GARCH UNIT ROOT TEST WITH 2 ENDOGENOUS BREAK DATES BY NARAYAN AND POPP 2010

I will be glad if you can assist me in running the code for this program.  i am using quarterly data from 2009 q1 to 2016q3.

kindly show me what i am to adjust from the codes to get the test result.  do not have a good background on gauss or coding but i have interest in running the code. thank you. below is the data with date i hope to use and the code also:

Code:

/* In return for using this codes, please cite the following relevant papers:/* In return for using this codes, please cite the following relevant papers:**** Paresh Kumar Narayan, Stephan Popp** A new Unit Root Test With Two Structural Breaks in Level and Slope at Unknown Time.** Journal of Applied Statistics, Volume 37, No.9, 1425-1428 (2010).****** Sequential procedure to estimate the break dates.**** Code written by Stephan Popp**** Copyright by Stephan Popp***/
new;FORMAT /M1 /ROS 8,4;library pgraph;_pdate="";_pcolor = 1;_pmcolor = ZEROS(8,1)|15;
/*path1 = "Z:\\Stephan.Popp\\Projekte\\perronmultiple\\gauss\\output.out";
OUTPUT FILE = ^path1 ON;*/
"******************************************";"date: " datestr(date) "  time:  " timestr(time);"*******************************************";
/* DATA**  the first column should contain the years or quarters (i.e. 1973.3)**  and the second the observations*/
load yy[231,1]=c:\data\roil.txt; @ change path settings and yy[here include # of observations ,2]@
/* Randomly chosen datasetyy=SEQA(1900,1,100)~RNDN(100,1);*/
XY(yy[.,1],yy[.,2]);
yyy= yy[.,2];kmax = 5; @ max lag @ttt = rows(yyy);tau = 0.2; @ trimming factor @
"Series: Log(.)";"Sample:" yy[1,1]~yy[rows(yyy),1];"# observations:" ttt;"maximum lag " kmax;"trimming factor " tau;
"program: popp2break.prg" ;"*******************************";tbunter = MAXC(3+kmax|CEIL(tau*ttt));             @ lower break date @"break date floor, effective tau " tbunter~yy[tbunter,1]~tbunter/ttt;tbober  = MINC(ttt-3-kmax|FLOOR((1-tau)*ttt));    @ upper break date @"break date ceiling, effective (1-tau) " tbober~yy[tbober,1]~tbober/ttt;"*******************************************";
/************** Model 0 *****************/
ergeb = ZEROS(2,1);     @ Vektor zum Speichern der Ergebnisse @tbopt = ZEROS(2,1);
ttb1 = tbunter;ttb2 = 0;do while ttb1 <= tbober; @ Schleife der durchlaufenden Bruchzeitpunkte T_B @     {rho1,trho1,ttheta1,pp1,rres1,varres1,tstat1,kk1} = mio2break0b(yyy,ttb1,ttb2,kmax);
if ABS(ttheta1) > ABS(ergeb[1,1]);   ergeb[1,1] = ttheta1;   ergeb[2,1] = ttb1; endif;
ttb1 = ttb1 + 1;endo;
tbopt[1] = ergeb[2,1];
ergeb = ZEROS(2,1);     @ Vektor zum Speichern der Ergebnisse @   ttb2 = tbunter;do while ttb2 <= tbober; if ABS(ttb2 - tbopt[1]) < 2; ttb2 = tbopt[1] + 2; endif;
{rho1,trho1,ttheta1,pp1,rres1,varres1,tstat1,kk1} = mio2break0b(yyy,tbopt[1],ttb2,kmax);
if ABS(ttheta1) > ABS(ergeb[1,1]);   ergeb[1,1] = ttheta1;   ergeb[2,1] = ttb2; endif;
ttb2 = ttb2 + 1;  endo;
tbopt[2] = ergeb[2,1];tbopt = SORTC(tbopt,1);
{rho1,trho1,ttheta1,pp1,rres1,varres1,tstat1,kk1} = mio2break0b(yyy,tbopt[1],tbopt[2],kmax);
"****************";"output";"****************";"model type M0";"first break  " tbopt[1]~yy[tbopt[1],1]~tbopt[1]/ROWS(yy);"second break " tbopt[2]~yy[tbopt[2],1]~tbopt[2]/ROWS(yy);"phi = rho-1  " rho1;"t value      " trho1;"optimal lag  " kk1;"variance     " varres1;"------";"coeff tstat  ";"yyverz constant du1verz du2verz dtb1 dtb2";pp1~tstat1;

/************** Model 1 *****************/
ergeb = ZEROS(2,1);     @ Vektor zum Speichern der Ergebnisse @tbopt = ZEROS(2,1);
ttb1 = tbunter;ttb2 = 0;do while ttb1 <= tbober; @ Schleife der durchlaufenden Bruchzeitpunkte T_B @     {rho1,trho1,ttheta1,pp1,rres1,varres1,tstat1,kk1} = mio2break1b(yyy,ttb1,ttb2,kmax);
if ABS(ttheta1) > ABS(ergeb[1,1]);   ergeb[1,1] = ttheta1;   ergeb[2,1] = ttb1; endif;
ttb1 = ttb1 + 1;endo;
tbopt[1] = ergeb[2,1];
ergeb = ZEROS(2,1);     @ Vektor zum Speichern der Ergebnisse @   ttb2 = tbunter;do while ttb2 <= tbober; if ABS(ttb2 - tbopt[1]) < 2; ttb2 = tbopt[1] + 2; endif;
{rho1,trho1,ttheta1,pp1,rres1,varres1,tstat1,kk1} = mio2break1b(yyy,tbopt[1],ttb2,kmax);
if ABS(ttheta1) > ABS(ergeb[1,1]);   ergeb[1,1] = ttheta1;   ergeb[2,1] = ttb2; endif;
ttb2 = ttb2 + 1;  endo;
tbopt[2] = ergeb[2,1];tbopt = SORTC(tbopt,1);
{rho1,trho1,ttheta1,pp1,rres1,varres1,tstat1,kk1} = mio2break1b(yyy,tbopt[1],tbopt[2],kmax);
"output";"****************";"model type M1";"first break  " tbopt[1]~yy[tbopt[1],1]~tbopt[1]/ROWS(yy);"second break " tbopt[2]~yy[tbopt[2],1]~tbopt[2]/ROWS(yy);"phi = rho-1  " rho1;"t value      " trho1;"optimal lag  " kk1;"variance     " varres1;"------";"coeff tstat  ";"yyverz constant du1verz du2verz dtb1 dtb2 time";pp1~tstat1;

/************** Model 2 *****************/
ergeb = ZEROS(2,1);     @ Vektor zum Speichern der Ergebnisse @tbopt = ZEROS(2,1);
ttb1 = tbunter;ttb2 = 0;do while ttb1 <= tbober; @ Schleife der durchlaufenden Bruchzeitpunkte T_B @     {rho1,trho1,ttheta1,pp1,rres1,varres1,tstat1,kk1} = mio2break2b(yyy,ttb1,ttb2,kmax);
if ABS(ttheta1) > ABS(ergeb[1,1]);   ergeb[1,1] = ttheta1;   ergeb[2,1] = ttb1; endif;
ttb1 = ttb1 + 1;endo;
tbopt[1] = ergeb[2,1];
ergeb = ZEROS(2,1);     @ Vektor zum Speichern der Ergebnisse @   ttb2 = tbunter;do while ttb2 <= tbober; if ABS(ttb2 - tbopt[1]) < 2; ttb2 = tbopt[1] + 2; endif;
{rho1,trho1,ttheta1,pp1,rres1,varres1,tstat1,kk1} = mio2break2b(yyy,tbopt[1],ttb2,kmax);
if ABS(ttheta1) > ABS(ergeb[1,1]);   ergeb[1,1] = ttheta1;   ergeb[2,1] = ttb2; endif;
ttb2 = ttb2 + 1;  endo;
tbopt[2] = ergeb[2,1];tbopt = SORTC(tbopt,1);
{rho1,trho1,ttheta1,pp1,rres1,varres1,tstat1,kk1} = mio2break2b(yyy,tbopt[1],tbopt[2],kmax);
"output";"****************";"model type M2";"first break  " tbopt[1]~yy[tbopt[1],1]~tbopt[1]/ROWS(yy);"second break " tbopt[2]~yy[tbopt[2],1]~tbopt[2]/ROWS(yy);"phi = rho-1  " rho1;"t value      " trho1;"optimal lag  " kk1;"variance     " varres1;"------";"coeff tstat  ";"yylagged constant du1lagged du2lagged dtb1 dtb2 dt1verz dt2verz time";pp1~tstat1;   OUTPUT OFF;
end;
proc(8)=mio2break0b(yy,ttb1,ttb2,kkmax);    local t,mu,du1,du1verz,du2,du2verz,dtb1,dtb2,yyverz,dyy,x1,x,pp,rres; local varres,varpp,tstat,kk,ampel,xx,ttheta;
t=ROWS(yy);    mu=ONES(t,1);              @ Absolutglied @ du1=ZEROS(t,1);             @ Dummy @ du1[ttb1+1:t,1]=ONES(t-ttb1,1); du1verz=lagn(du1,1); @ Dummy verz?gert DU_t-1 @ dtb1=ZEROS(t,1);            @ Impulsdummy @ dtb1[ttb1+1:ttb1+1,1]=1; yyverz=lagn(yy,1); dyy=yy-yyverz;
if ttb2 == 0;   x1=yyverz~mu~du1verz~dtb1; else;   du2=ZEROS(t,1);             @ Dummy @   du2[ttb2+1:t,1]=ONES(t-ttb2,1);   du2verz=lagn(du2,1); @ Dummy verz?gert DU_t-1 @   dtb2=ZEROS(t,1);            @ Impulsdummy @   dtb2[ttb2+1:ttb2+1,1]=1;
x1=yyverz~mu~du1verz~du2verz~dtb1~dtb2; endif;   if kkmax == 0;   x = TRIMR(x1,1,0);   dyy = TRIMR(dyy,1,0);    pp = INV(x'x)*x'dyy; @ Sch?tzvektor p @   rres = dyy - x*pp; @ Residuenvektor @      varres =(rres'*rres)/(ROWS(x)-COLS(x));                 @ Sch?tzung der Residualvarianz @      varpp = DIAG(rres'rres*INV(x'x)/(ROWS(x)-COLS(x))); @ Sch?tzung der Varianz-Kovarianzmatrix @      tstat = pp./SQRT(varpp);                                  @ t-Statistiken der Parametersch?tzungen @   kk = 0; else;   x = TRIMR(x1~SHIFTR(dyy',SEQA(1,1,kkmax),0)',kkmax+1,0);   dyy = TRIMR(dyy,kkmax+1,0);   kk = kkmax;   ampel = 0;
do while kk >= 0 and ampel == 0;     xx = x[.,1:COLS(x1)+kk];
pp = INV(xx'xx)*xx'dyy; @ Sch?tzvektor p @     rres = dyy - xx*pp; @ Residuenvektor @        varres =(rres'*rres)/(ROWS(xx)-COLS(xx));                 @ Sch?tzung der Residualvarianz @        varpp = DIAG(rres'rres*INV(xx'xx)/(ROWS(xx)-COLS(xx))); @ Sch?tzung der Varianz-Kovarianzmatrix @        tstat = pp./SQRT(varpp);                                  @ t-Statistiken der Parametersch?tzungen @       if ABS(tstat[COLS(x1)+kk]) > 1.96; ampel = 1; endif;   @ Signifkanzniveau 10%: 1.96 @     kk = kk-1;   endo;   kk = kk+1;
endif;
if ttb2 == 0;   ttheta = tstat[4]; else;   ttheta = tstat[6]; endif;   retp(pp[1],tstat[1],ttheta,pp,rres,varres,tstat,kk);
endp;
proc(8)=mio2break1b(yy,ttb1,ttb2,kkmax);    local t,mu,zeit,du1,du1verz,du2,du2verz,dtb1,dtb2,yyverz,dyy,x1,x,pp,rres; local varres,varpp,tstat,kk,ampel,xx,ttheta;
t=ROWS(yy);    mu=ONES(t,1);              @ Absolutglied @    zeit=SEQA(1,1,t);          @ Zeittrend @ du1=ZEROS(t,1);             @ Dummy @ du1[ttb1+1:t,1]=ONES(t-ttb1,1); du1verz=lagn(du1,1); @ Dummy verz?gert DU_t-1 @ dtb1=ZEROS(t,1);            @ Impulsdummy @ dtb1[ttb1+1:ttb1+1,1]=1; yyverz=lagn(yy,1); dyy=yy-yyverz;
if ttb2 == 0;   x1=yyverz~mu~du1verz~dtb1~zeit; else;   du2=ZEROS(t,1);             @ Dummy @   du2[ttb2+1:t,1]=ONES(t-ttb2,1);   du2verz=lagn(du2,1); @ Dummy verz?gert DU_t-1 @   dtb2=ZEROS(t,1);            @ Impulsdummy @   dtb2[ttb2+1:ttb2+1,1]=1;
x1=yyverz~mu~du1verz~du2verz~dtb1~dtb2~zeit; endif;   if kkmax == 0;   x = TRIMR(x1,1,0);   dyy = TRIMR(dyy,1,0);    pp = INV(x'x)*x'dyy; @ Sch?tzvektor p @   rres = dyy - x*pp; @ Residuenvektor @      varres =(rres'*rres)/(ROWS(x)-COLS(x));                 @ Sch?tzung der Residualvarianz @      varpp = DIAG(rres'rres*INV(x'x)/(ROWS(x)-COLS(x))); @ Sch?tzung der Varianz-Kovarianzmatrix @      tstat = pp./SQRT(varpp);                                  @ t-Statistiken der Parametersch?tzungen @   kk = 0; else;   x = TRIMR(x1~SHIFTR(dyy',SEQA(1,1,kkmax),0)',kkmax+1,0);   dyy = TRIMR(dyy,kkmax+1,0);   kk = kkmax;   ampel = 0;
do while kk >= 0 and ampel == 0;     xx = x[.,1:COLS(x1)+kk];
pp = INV(xx'xx)*xx'dyy; @ Sch?tzvektor p @     rres = dyy - xx*pp; @ Residuenvektor @        varres =(rres'*rres)/(ROWS(xx)-COLS(xx));                 @ Sch?tzung der Residualvarianz @        varpp = DIAG(rres'rres*INV(xx'xx)/(ROWS(xx)-COLS(xx))); @ Sch?tzung der Varianz-Kovarianzmatrix @        tstat = pp./SQRT(varpp);                                  @ t-Statistiken der Parametersch?tzungen @       if ABS(tstat[COLS(x1)+kk]) > 1.96; ampel = 1; endif;   @ Signifkanzniveau 10%: 1.96 @     kk = kk-1;   endo;   kk = kk+1;
endif;
if ttb2 == 0;   ttheta = tstat[4]; else;   ttheta = tstat[6]; endif;   retp(pp[1],tstat[1],ttheta,pp,rres,varres,tstat,kk);
endp;
proc(8)=mio2break2b(yy,ttb1,ttb2,kkmax);    local t,mu,zeit,du1,du1verz,du2,du2verz,dtb1,dtb2,dt1,dt1verz,dt2,dt2verz,yyverz,dyy,x1,x,pp,rres; local varres,varpp,tstat,kk,ampel,xx,ttheta;
t=ROWS(yy);    mu=ONES(t,1);              @ Absolutglied @    zeit=SEQA(1,1,t);          @ Zeittrend @ du1=ZEROS(t,1);             @ Dummy @ du1[ttb1+1:t,1]=ONES(t-ttb1,1); du1verz=lagn(du1,1); @ Dummy verz?gert DU_t-1 @ dtb1=ZEROS(t,1);            @ Impulsdummy @ dtb1[ttb1+1:ttb1+1,1]=1; dt1 = ZEROS(t,1); dt1[ttb1+1:t,1] = SEQA(1,1,(t-ttb1)); dt1verz=lagn(dt1,1); yyverz=lagn(yy,1); dyy=yy-yyverz;
if ttb2 == 0;   x1=yyverz~mu~du1verz~dtb1~dt1verz~zeit; else;   du2=ZEROS(t,1);             @ Dummy @   du2[ttb2+1:t,1]=ONES(t-ttb2,1);   du2verz=lagn(du2,1); @ Dummy verz?gert DU_t-1 @   dtb2=ZEROS(t,1);            @ Impulsdummy @   dtb2[ttb2+1:ttb2+1,1]=1;   dt2 = ZEROS(t,1);   dt2[ttb2+1:t,1] = SEQA(1,1,(t-ttb2));   dt2verz=lagn(dt2,1);
x1=yyverz~mu~du1verz~du2verz~dtb1~dtb2~dt1verz~dt2verz~zeit; endif;   if kkmax == 0;   x = TRIMR(x1,1,0);   dyy = TRIMR(dyy,1,0);    pp = INV(x'x)*x'dyy; @ Sch?tzvektor p @   rres = dyy - x*pp; @ Residuenvektor @      varres =(rres'*rres)/(ROWS(x)-COLS(x));                 @ Sch?tzung der Residualvarianz @      varpp = DIAG(rres'rres*INV(x'x)/(ROWS(x)-COLS(x))); @ Sch?tzung der Varianz-Kovarianzmatrix @      tstat = pp./SQRT(varpp);                                  @ t-Statistiken der Parametersch?tzungen @   kk = 0; else;   x = TRIMR(x1~SHIFTR(dyy',SEQA(1,1,kkmax),0)',kkmax+1,0);   dyy = TRIMR(dyy,kkmax+1,0);   kk = kkmax;   ampel = 0;
do while kk >= 0 and ampel == 0;     xx = x[.,1:COLS(x1)+kk];
pp = INV(xx'xx)*xx'dyy; @ Sch?tzvektor p @     rres = dyy - xx*pp; @ Residuenvektor @        varres =(rres'*rres)/(ROWS(xx)-COLS(xx));                 @ Sch?tzung der Residualvarianz @        varpp = DIAG(rres'rres*INV(xx'xx)/(ROWS(xx)-COLS(xx))); @ Sch?tzung der Varianz-Kovarianzmatrix @        tstat = pp./SQRT(varpp);                                  @ t-Statistiken der Parametersch?tzungen @       if ABS(tstat[COLS(x1)+kk]) > 1.96; ampel = 1; endif;   @ Signifkanzniveau 10%: 1.96 @     kk = kk-1;   endo;   kk = kk+1;
endif;
if ttb2 == 0;   ttheta = tstat[4]; else;   ttheta = tstat[6]; endif;   retp(pp[1],tstat[1],ttheta,pp,rres,varres,tstat,kk);
endp;
proc diff(x,k) ;     if ( k == 0) ;        retp(x) ;     endif ;     retp(trimr(x,k,0)-trimr(lagn(x,k),k,0)) ;endp ;
proc lagn(x,n);    local y;    y = shiftr(x', n, (miss(0, 0))');    retp(y');endp;