==AUTOGARCH Autogarch setup template /$ /$ Illustrates Autogarch Setup /$ User supplied series /$ b34sexec options ginclude('gas.b34'); b34srun; b34sexec matrix; call loaddata; /$ /$ Settings here /$ series=gasout; doplot =0; doplot2=0; nacf=24; npacf=24; s=2.0; seas=0; /$ for forecasts need to change here. Default = 5 at end /$ For "custom" setupsturn on indicated commands nf =5; nt =norows(series); nf2=5; nt2=norows(series); /$ /$ rtest has residual and y plots /$ call load(rtest); call load(rtest2); /$ /$ This roottol setting forces no differencing /$ /$ call autobj(gasout :print :nac 24 :npac 24 /$ :roottol .99 :autobuild ); /$ This turns off differencing call autobj(series :print :nac nacf :npac npacf /$ :nodif /$ :printsteps :seasonal seas /$ :rawacfpacf /$ :difrawacf /$ :rdif /$ :sdif :forecast index(nf nt) :spiketol s :autobuild ); if(doplot.ne.0) call rtest(%res,series,nacf); if(doplot2.ne.0)call rtest2(%res,series,nacf); /$ Default let program decide ressq=%res*%res; call autobj(ressq :print :nac nacf :npac npacf /$ :printsteps :nodif /$ :rawacfpacf /$ :difrawacf :forecast index(nf2 nt2) :spiketol s :autobuild ); if(doplot.ne.0) call rtest2(%res,rrsq,48); if(doplot2.ne.0)call rtest2(%res,rrsq,48); b34srun; == ==BOOST Boost & Boost2 tests /; /; Implements traditional "centered" boosting. For non-centered /; OLS boosting see BOOST3 and BOOST4 for forecasting /; b34sexec options ginclude('b34sdata.mac') member(efron_1); b34srun; /; b34sexec options ginclude('b34sdata.mac') member(efron_2); b34srun; /; /; Implements OLS boosting as discribed in "Least Angle Regression" /; by Bradley Efron, Trevor Hastie, Iain Johnson and Robert Tibshirani /; The Annals of Statistics Vol 32 No. 2 (April 2004) pp 407-451 /; /; Modified boosting also available if iboost set = 1 /; b34sexec matrix; call loaddata; call load(center :staging); call load(center2 :staging); call load(boost :staging); call echooff; /; iboost=0 => ols boost /; iboost=1 => modified OLS boost iboost=0; /; iboost=1; e=.5; ntry=100; /; itype=0 => ols /; itype=1 => mars /; itype=2 => gam /; itype=3 => l1 /; itype=4 => minimax itype=0; iprint=1; x=center2(catcol(age sex bmi bp s1 s2 s3 s4 s5 s6)); y=y-mean(y); /; Set up the base case. if(itype.eq.0)call olsq(y,x :noint :print); if(itype.eq.1)call marspline(y x :print :nk 20); if(itype.eq.2)call gamfit(y x[predictor,3] :print); if(itype.eq.3)call olsq(y,x :noint :print :l1); if(itype.eq.4)call olsq(y,x :noint :print :minimax); call print(ccf(%y,%yhat)); base1=ccf(%y,%yhat); fit=array(ntry:); do i=1,ntry; if(iboost.eq.0)call boost( y,yhat,res,x, in,e,i,itype,iprint); if(iboost.eq.1)call boost2(y,yhat,res,x,xbuild,in,e,i,itype,iprint); /; if(iprint.ne.0)call graph(y,yhat); fit(i)=ccf(y,yhat); call outstring(1,2,'Iteration'); call outinteger(20,2,i); call outstring(1,3,'Correlation'); call outdouble(20,3,fit(i)); enddo; base=array(norows(fit):)+base1; if(itype.eq.0)call char1(cc1,'OLS' ); if(itype.eq.1)call char1(cc1,'MARS' ); if(itype.eq.2)call char1(cc1,'GAM' ); if(itype.eq.3)call char1(cc1,'L1' ); if(itype.eq.4)call char1(cc1,'MINIMAX'); cc=' Boosting based Correlation of y and yhat given eps ='; if(iboost.eq.1) cc=' Modified Boosting based Correlation of y and yhat given eps ='; call ir8tostr(e,cc2,'(f8.4)'); cc=catrow(cc1,cc,cc2); call graph(base,fit :heading cc :file 'boost.wmf' :nocontact :nolabel :pgborder); call print(fit); b34srun; == ==BOOST3 OLS BOOSTING Model - HHS Mods with Forecasting b34sexec options ginclude('b34sdata.mac') member(efron_1); b34srun; /; b34sexec options ginclude('b34sdata.mac') member(efron_2); b34srun; /; /; Implements OLS boosting as discribed in "Least Angle Regression" /; by Bradley Efron, Trevor Hastie, Iain Johnson and Robert Tibshirani /; The Annals of Statistics Vol 32 No. 2 (April 2004) pp 407-451 /; b34sexec matrix; call loaddata; call load(center :staging); call load(center2 :staging); call load(boost :staging); call echooff; /; itype=0 => ols /; itype=1 => mars not ready /; itype=2 => gam not ready /; itype=3 => L1 /; itype=4 => MINIMAX itype=0; iprint=0; e=.5 ; ntry=1000; x=catcol(age sex bmi bp s1 s2 s3 s4 s5 s6); /; Gets the base if(itype.eq.0)call olsq(y,x :print); if(itype.eq.1)call marspline(y x :print :nk 20); if(itype.eq.2)call gamfit(y x[predictor,3] :print); if(itype.eq.3)call olsq(y,x :print :l1); if(itype.eq.4)call olsq(y,x :print :minimax); call print(ccf(%y,%yhat)); base1=ccf(%y,%yhat); fit =array(ntry:); beta1 =array(ntry:); beta2 =array(ntry:); in =idint(array(ntry:)); do i=1,ntry; call boost3(y,yhat,res,x,in,beta1,beta2,e,i,itype,iprint); fit(i)=ccf(y,yhat); call outstring(1,2,'Iteration'); call outinteger(20,2,i); call outstring(1,3,'Correlation'); call outstring(1,3,'Correlation'); call outdouble(20,3,fit(i)); enddo; call tabulate(in,fit,beta1,beta2); /; Testing if can "forecast" using saved Model call boost4(yhat2,x,in,beta1,beta2,e); call tabulate(yhat,yhat2); base=array(norows(fit):)+base1; if(itype.eq.0)call char1(cc1,'OLS' ); if(itype.eq.1)call char1(cc1,'MARS' ); if(itype.eq.2)call char1(cc1,'GAM' ); if(itype.eq.3)call char1(cc1,'L1' ); if(itype.eq.4)call char1(cc1,'MINIMAX'); cc=' HHS Boosting based Correlation of y and yhat given eps ='; call ir8tostr(e,cc2,'(f8.4)'); cc=catrow(cc1,cc,cc2); call graph(base,fit :heading cc :file 'boost.wmf' :nocontact :nolabel :pgborder); call print(fit); b34srun; == ==BOOST4 OLS BOOSTING Model - HHS Mods with Forecasting b34sexec options ginclude('b34sdata.mac') member(efron_1); b34srun; /; b34sexec options ginclude('b34sdata.mac') member(efron_2); b34srun; /; /; Implements OLS boosting as discribed in "Least Angle Regression" /; by Bradley Efron, Trevor Hastie, Iain Johnson and Robert Tibshirani /; The Annals of Statistics Vol 32 No. 2 (April 2004) pp 407-451 /; b34sexec matrix; call loaddata; call load(center :staging); call load(center2 :staging); call load(boost :staging); call echooff; /; itype=0 => ols /; itype=1 => mars not ready /; itype=2 => gam not ready /; itype=3 => L1 /; itype=4 => MINIMAX itype=0; iprint=0; e=.5 ; ntry=1000; x=catcol(age sex bmi bp s1 s2 s3 s4 s5 s6); /; Gets the base if(itype.eq.0)call olsq(y,x :print); if(itype.eq.1)call marspline(y x :print :nk 20); if(itype.eq.2)call gamfit(y x[predictor,3] :print); if(itype.eq.3)call olsq(y,x :print :l1); if(itype.eq.4)call olsq(y,x :print :minimax); call print(ccf(%y,%yhat)); base1=ccf(%y,%yhat); fit =array(ntry:); beta1 =array(ntry:); beta2 =array(ntry:); in =idint(array(ntry:)); do i=1,ntry; call boost3(y,yhat,res,x,in,beta1,beta2,e,i,itype,iprint); fit(i)=ccf(y,yhat); call outstring(1,2,'Iteration'); call outinteger(20,2,i); call outstring(1,3,'Correlation'); call outstring(1,3,'Correlation'); call outdouble(20,3,fit(i)); enddo; call tabulate(in,fit,beta1,beta2); /; Testing if can "forecast" using saved Model call boost4(yhat2,x,in,beta1,beta2,e); call tabulate(yhat,yhat2); base=array(norows(fit):)+base1; if(itype.eq.0)call char1(cc1,'OLS' ); if(itype.eq.1)call char1(cc1,'MARS' ); if(itype.eq.2)call char1(cc1,'GAM' ); if(itype.eq.3)call char1(cc1,'L1' ); if(itype.eq.4)call char1(cc1,'MINIMAX'); cc=' HHS Boosting based Correlation of y and yhat given eps ='; call ir8tostr(e,cc2,'(f8.4)'); cc=catrow(cc1,cc,cc2); call graph(base,fit :heading cc :file 'boost.wmf' :nocontact :nolabel :pgborder); call print(fit); b34srun; == ==CATNAME Build an augmented NAME vector /; /; Application of LTS to MARS and GAM Models /; b34sexec options ginclude('gas.b34'); b34srun; b34sexec matrix; call loaddata; call load(lts :staging); call load(upmatrix :staging); call load(catname :staging); call echooff; n=6; p=.9; iprint=1; call olsq(gasout gasin{1 to n} gasout{1 to n} :print :savex); call lts(%y,%x,%res,%names,%lag,%coef,%se,%t,p,newy,newx,ihold, iprint); call marspline(newy newx :print :mi 2 :nk 20 :holdout ihold); call upmatrix(newx,%names,%lag,level(),newnames); add ='[predictor,3]'; /; /; Do not update constant /; call catname(newnames,add,augname,1); call gamfit(newy argument(augname) :holdout ihold :print); b34srun; == ==CATNAME2 Example of just a name update b34sexec matrix; call loaddata; call load(catname :staging); call echooff; testn=' houston'; call character(add,'.dat_99999999'); call print(testn,add); call catname(testn,add,rhs,0); call print(rhs); b34srun; == ==CATNAME3 Various examples of CATNAME b34sexec matrix; call load(catname :staging); call echooff; call character(t1, '6 thisis'); call character(t2 '.dat_123'); t2='.dat'; call catname(t1,t2,rhs,0); call print(t1,t2); call print(rhs); call names(all); testn=' houston_is_xxx'; call catname(testn,t2,new2,0); call print( ' houston_is_xxx.dat', t2,testn,new2); v=namelist(don,sue,joe); call character(t2,'.dat_12345678'); call print(testn,t2); /; call catname(testn,t2,rhs,0); call print(rhs); call catname(v,t2,rhs,0); call print(rhs); call catname(v,t2,rhs,1); call print(rhs); b34srun; == ==CENTER Center a series or 2 D object b34sexec options ginclude('b34sdata.mac') member(efron_1); b34srun; b34sexec matrix; call loaddata; call load(center :staging); call echooff; y=center(y); call print('Test - Mean of centered variable ', mean(y)); call print('Test - Variance of centered variable ',variance(y)); xx=catcol(age,sex,bmi,bp,s1,s2,s3,s4,s5,s6); cxx=center(xx); mm=array(10:); vv=array(10:); do i=1,10; mm(i)=mean( cxx(,i)); vv(i)=variance(cxx(,i)); enddo; call print('Test = Mean of centered matrix ', mm); call print('Test = variance of centered matrix ', vv); b34srun; == ==CENTER2 Standardize a Series mean=0. Unit length. b34sexec options ginclude('b34sdata.mac') member(efron_1); b34srun; b34sexec matrix; call loaddata; call load(center :staging); call load(center2:staging); call echooff; call print('Centering such that Mean=0. Variance = 1':); x=center(catcol(age sex bmi bp s1 s2 s3 s4 s5 s6)); call print(mean(x(,1)),variance(x(,1))); call print(mean(x(,2)),variance(x(,2))); call print(mean(x(,3)),variance(x(,3))); call print(mean(x(,4)),variance(x(,4))); call print(mean(x(,5)),variance(x(,5))); call print(mean(x(,6)),variance(x(,6))); call print(mean(x(,7)),variance(x(,7))); call print(mean(x(,8)),variance(x(,8))); call print(mean(x(,9)),variance(x(,9))); call print(mean(x(,10)),variance(x(,10))); call print(' ':); call print('Standardizing such that Mean=0 sumsq(series)=1':); x=center2(catcol(age sex bmi bp s1 s2 s3 s4 s5 s6)); call print(mean(x(,1)),variance(x(,1))); call print(mean(x(,2)),variance(x(,2))); call print(mean(x(,3)),variance(x(,3))); call print(mean(x(,4)),variance(x(,4))); call print(mean(x(,5)),variance(x(,5))); call print(mean(x(,6)),variance(x(,6))); call print(mean(x(,7)),variance(x(,7))); call print(mean(x(,8)),variance(x(,8))); call print(mean(x(,9)),variance(x(,9))); call print(mean(x(,10)),variance(x(,10))); x1=x(,1); call print('Testing x1 ',sumsq(x1)); b34srun; b34sexec options ginclude('b34sdata.mac') member(efron_2); b34srun; == ==CGETR16 Easy way to read real*16 data into matrix /; /; Easy way to read data into real*16 or VPA /; b34sexec matrix; datacards; 10.07E0 77.6E0 14.73E0 114.9E0 17.94E0 141.1E0 23.93E0 190.8E0 29.61E0 239.9E0 35.18E0 289.0E0 40.02E0 332.8E0 44.82E0 378.4E0 50.76E0 434.8E0 55.05E0 477.3E0 61.01E0 536.8E0 66.40E0 593.1E0 75.47E0 689.1E0 81.78E0 760.0E0 b34sreturn; call echooff; call load(getvpa :staging); call load(cgetr16 :staging); call load(cgetvpa :staging); call print(' ':); call print('Test cases with 100% real*16 and VPA data loading':); call print('_________________________________________________':); call cgetr16(14,2,xx); call print(xx); call cgetvpa(14,2,xvpa); call print(xvpa); b34srun; == ==CGETVPA Easy way to read VPA data into Matrix /; /; Easy way to read data into real*16 or VPA /; b34sexec matrix; datacards; 10.07E0 77.6E0 14.73E0 114.9E0 17.94E0 141.1E0 23.93E0 190.8E0 29.61E0 239.9E0 35.18E0 289.0E0 40.02E0 332.8E0 44.82E0 378.4E0 50.76E0 434.8E0 55.05E0 477.3E0 61.01E0 536.8E0 66.40E0 593.1E0 75.47E0 689.1E0 81.78E0 760.0E0 b34sreturn; call echooff; call load(getvpa :staging); call load(cgetr16 :staging); call load(cgetvpa :staging); call print(' ':); call print('Test cases with 100% real*16 and VPA data loading':); call print('_________________________________________________':); call cgetr16(14,2,xx); call print(xx); call cgetvpa(14,2,xvpa); call print(xvpa); b34srun; == ==COR Correlation Test Cases b34sexec options ginclude('gas.b34'); b34srun; b34sexec matrix; call loaddata; call load(cov :staging); call load(cor :staging); call echooff; h=catcol(gasin,gasout); /; call print(h); call print('Using Built in Command for Correlation':); call print(ccf(h)); call print('Covariance using function cov':); call print(cov(h)); call print('Correlation using function cor':); call print(cor(h)); h16=r8tor16(h); call print('Covariance using function cov':); call print(cov(h16)); call print('Correlation using function cor':); call print(cor(h16)); h_vpa=vpa(h); call print('Covariance using function cov':); call print(cov(h_vpa)); call print('Correlation using function cor':); call print(cor(h_vpa)); b34srun; == ==COND Matrix Condition /; Illustrates norm and cond /; Optionally Matlab can be called /; B34S tracks 100% Matlab / LAPACK /; %b34slet domatlab=0; b34sexec options ginclude('gas.b34'); b34srun; b34sexec matrix; call loaddata; call load(cov :staging); call load(cor :staging); call load(norm :staging); call load(cond :staging); call echooff; lag=3; /; call olsq(gasout gasout{1 to lag} gasin{1 to lag} :print :savex); call print('Tests on Real*8 Cholesky':); call print('++++++++++++++++++++++++':); xpx=transpose(%x)*%x; call print('Norm 1 xpx ',norm(xpx,1):); call print('Norm 2 xpx ',norm(xpx,2):); call print('Norm 3 xpx ',norm(xpx,3):); call print('Norm 4 xpx ',norm(xpx,4):); call print('Cond 1 xpx ',cond(xpx,1):); call print('Cond 2 xpx ',cond(xpx,2):); call print('Cond 3 xpx ',cond(xpx,3):); call print('Cond 4 xpx ',cond(xpx,4):); call print('Norm 1 x ',norm(%x,1):); call print('Norm 2 x ',norm(%x,2):); call print('Norm 3 x ',norm(%x,3):); call print('Norm 4 x ',norm(%x,4):); call print('Cond 2 x ',cond(%x,2):); test=inv(xpx,testcond :gmat); testcond=1.0/testcond; call print('Condition from LAPACK ',testcond:); test=inv(xpx,testcond ); testcond=1.0/testcond; call print('Condition from LINPACK ',testcond:); %b34sif(&domatlab.ne.0)%then; call makematlab(%x :file 'xdata.m'); call load(rmatlab); call rmatlab; pgmcards; format long g x=getb34s('xdata.m'); xpx=x'*x; % ' disp('norm(xpx,1)') disp( norm(xpx,1) ) disp('norm(xpx,2)') disp( norm(xpx,2) ) disp('norm(xpx,inf)') disp( norm(xpx,inf) ) disp('norm(xpx,"fro")') disp( norm(xpx,'fro')) disp('cond(xpx,1)') disp( cond(xpx,1) ) disp('cond(xpx,2)') disp( cond(xpx,2) ) disp('cond(xpx,inf)') disp( cond(xpx,inf) ) disp('cond(xpx,"fro")') disp( cond(xpx,'fro') ) disp('norm(x,1)') disp( norm(x,1) ) disp('norm(x,2)') disp( norm(x,2) ) disp('norm(x,inf)') disp( norm(x,inf) ) disp('norm(x,"fro")') disp( norm(x,'fro') ) quit b34sreturn; %b34sendif; b34srun; == ==COR_2 Effect of precision on OLS Assumptions b34sexec options ginclude('gas.b34'); b34srun; b34sexec matrix; call loaddata; call load(cov :staging); call load(cor :staging); call echooff; lag=3; /; /; This options takes space /; dovpa=0; doreal4=1; call olsq(gasout gasout{1 to lag} gasin{1 to lag} :print :savex :diag); call print('Tests on Real*8 Cholesky':); call print('++++++++++++++++++++++++':); i=nocols(%x); test=%x; test(1,i)=%res; call print('Last Column is the residual':); call print(cor(test)); call olsq(gasout gasout{1 to lag} gasin{1 to lag} :print :savex :qr :diag); call print('Tests on Real*8 QR':); call print('++++++++++++++++++':); i=nocols(%x); test=%x; test(1,i)=%res; call print('Last Column is the residual':); call print(cor(test)); %yhold=%y; %xhold=%x; %ytest=r8tor16(%yhold); %xtest=r8tor16(%xhold); call olsq(%ytest %xtest :print :diag :noint :savex); call print('Tests on Real*16 Cholesky':); call print('+++++++++++++++++++++++++':); i=nocols(%x); test=%x; test(1,i)=%res; call print('Last Column is the residual':); call print(cor(test)); call olsq(%ytest %xtest :print :diag :noint :qr :savex); call print('Tests on Real*16 QR':); i=nocols(%x); test=%x; test(1,i)=%res; call print('Last Column is the residual':); call print(cor(test)); if(dovpa.ne.0)then; call print('VPA Tests ':); call print('_________ ':); %ytest=vpa(%yhold); %xtest=vpa(%xhold); beta = inv(transpose(%xtest)*%xtest)*transpose(%xtest)*%ytest; call print('beta from VPA ':); call print(beta); i=nocols(%xtest); %res=%xtest*beta-%ytest; test=%xtest; test(1,i)=%res; call print('Last Column is the residual':); call print(cor(test)); endif; /; This example show how key method is when data is not /; saved with enough precision if(doreal4.ne.0)then; call print('Real*4 Tests ':); call print('_________ ':); %ytest=r8tor4(%yhold); %xtest=r8tor4(%xhold); beta = inv(transpose(%xtest)*%xtest)*transpose(%xtest)*%ytest; call print('beta from Real*4 using inverse ':); call print(beta); i=nocols(%xtest); %res=%xtest*beta-%ytest; test=%xtest; test(1,i)=%res; call print('Last Column is the residual':); call print(cor(test)); /; QR approach r4_r=qr( %xtest ,r4_q); /; r8_r=qr(r4tor8(%xtest),r8_q); /; call print('R looked at from real 4 to real 8 and from real 4', /; r8_r,r4_r); r4_yhat=r4_q*transpose(r4_q)*%ytest; r4_res =%ytest-r4_yhat; /; /; Beta not needed to get yhat from QR !! /; r4_beta=inv(r4_r)*transpose(r4_q)*%ytest; call print('beta from Real*4 using QR ':); call print(r4_beta); i=nocols(%xtest); test=%xtest; test(1,i)=r4_res; call print('Last Column is the residual':); call print(cor(test)); endif; b34srun; == ==COR_3 Effect of Data Precision on Accuracy /; /; Charles Renfro Test Regression /; b34sexec options ginclude('b34sdata.mac') member(grunfeld_4); b34srun; b34sexec matrix; call loaddata; call echooff; /; /; Set parameters for tests /; /; call print(' ':); call print('GE Equation':); call print('-----------':); call print(' ':); call print(' ':); call load(cov :staging); call load(cor :staging); call echooff; /; dovpa=1; call olsq(gei gef gec :print :savex); call print('Tests on Real*8 Cholesky':); i=nocols(%x); test=%x; test(1,i)=%res; call print('Last Column is the residual':); call print(cor(test)); call olsq(gei gef gec :print :savex :qr); call print('Tests on Real*8 QR':); i=nocols(%x); test=%x; test(1,i)=%res; call print('Last Column is the residual':); call print(cor(test)); %yhold=%y; %xhold=%x; %ytest=r8tor16(%yhold); %xtest=r8tor16(%xhold); call olsq(%ytest %xtest :print :noint :savex); call print('Tests on Real*16 Cholesky':); i=nocols(%x); test=%x; test(1,i)=%res; call print('Last Column is the residual':); call print(cor(test)); call olsq(%ytest %xtest :print :noint :qr :savex); call print('Tests on Real*16 QR':); i=nocols(%x); test=%x; test(1,i)=%res; call print('Last Column is the residual':); call print(cor(test)); /; /; This example show how key method is when data is not /; saved with enough precision call print('Real*4 Tests ':); call print('_________ ':); %ytest=r8tor4(%yhold); %xtest=r8tor4(%xhold); beta = inv(transpose(%xtest)*%xtest)*transpose(%xtest)*%ytest; call print('beta from Real*4 using inverse ':); call print(beta); i=nocols(%xtest); %res=%xtest*beta-%ytest; test=%xtest; test(1,i)=%res; call print('Last Column is the residual':); call print(cor(test)); /; QR approach r4_r=qr( %xtest ,r4_q); /; r8_r=qr(r4tor8(%xtest),r8_q); /; call print('R looked at from real 4 to real 8 and from real 4', /; r8_r,r4_r); r4_yhat=r4_q*transpose(r4_q)*%ytest; r4_res =%ytest-r4_yhat; /; /; Beta not needed to get yhat from QR !! /; r4_beta=inv(r4_r)*transpose(r4_q)*%ytest; call print('beta from Real*4 using QR ':); call print(r4_beta); i=nocols(%xtest); test=%xtest; test(1,i)=r4_res; call print('Last Column is the residual':); call print(cor(test)); if(dovpa.ne.0)then; call print('VPA Tests ':); call print('_________ ':); %ytest=vpa(%yhold); %xtest=vpa(%xhold); beta = inv(transpose(%xtest)*%xtest)*transpose(%xtest)*%ytest; call print('beta from VPA ':); call print(beta); i=nocols(%xtest); %res=%xtest*beta-%ytest; test=%xtest; test(1,i)=%res; call print('Last Column is the residual':); call print(cor(test)); endif; b34srun; == ==COV Covariance/Correlation Test Cases b34sexec options ginclude('gas.b34'); b34srun; b34sexec matrix; call loaddata; call load(cov :staging); call load(cor :staging); call echooff; h=catcol(gasin,gasout); /; call print(h); call print('Covariance using function cov':); call print(cov(h)); call print('Correlation using function cor':); call print(cor(h)); call print('Real*16 results':); h16=r8tor16(h); call print('Covariance using function cov':); call print(cov(h16)); call print('Correlation using function cor':); call print(cor(h16)); b34srun; == ==COV2 Tests of cov2( ) /; /; This illustrates two ways to get the variance-Covariance Matrix /; Using Greene's Idempotentent Matrix and real*4. If the data is not /; scaled, there can be problems of accuracy that are detected using /; real*16 results as the benchmark /; b34sexec options ginclude('gas.b34'); b34srun; b34sexec matrix; call loaddata; call load(cov :staging); call load(cov2 :staging); call load(cor :staging); call echooff; scale=1000.; h=catcol(gasin,scale*gasout); /; call print(h); call print('Covariance using function cov':); call print(cov(h)); call print('Covariance using function cov2':); call print(cov2(h)); call print('Difference of two methods using real*8':); call print(cov(h)-cov2(h)); call print('Correlation using function cor':); call print(cor(h)); call print('Real*16 results':); h16=r8tor16(h); call print('Covariance using function cov':); call print(cov(h16)); call print('Covariance using function cov2':); call print(cov2(h16)); call print('Difference of the two methods using real*16 ':); call print(cov(h)-cov2(h)); /; Testing which is closer real8_1=afam(cov(h)); real8_2=afam(cov2(h)); call print('Difference against real16 for real8_1 & real8_2 ':); call print('Where Traditional Method = real8_1 ':); call print('Where M0 Method = real8_2 ':); call print(r8tor16(real8_1)-afam(cov(h16))); call print(r8tor16(real8_2)-afam(cov(h16))); call print('Correlation using function cor':); call print(cor(h16)); /; real*4 results h4=r8tor4(h); call print(' ':); call print('Where Traditional Method = real4_1 ':); call print('Where M0 Method = real4_2 ':); real4_1=afam(cov(h4)); real4_2=afam(cov2(h4)); call print(real4_1,real4_2); call print('Difference against real16 for real4_1 & real4_2 ':); call print(r8tor16(r4tor8(real4_1))-afam(cov(h16))); call print(r8tor16(r4tor8(real4_2))-afam(cov(h16))); call print('Correlation using function cor':); call print(cor(h4)); b34srun; == ==COVTEST Tests against Variance-Covarience Matrix b34sexec options ginclude('gas.b34'); b34srun; b34sexec data set cov; b34srun; b34sexec matrix; call echooff; call loaddata; call load(cov :staging); call load(cov2 :staging); z=catcol(time,gasin gasout,constant); call print('Traditional Method ',cov(z)); call print('Idempotent Matrix Method',cov2(z)); b34srun; == ==DDOT Vector Product with an offset b34sexec options ginclude('greene.mac') member(nf5_1); b34srun; b34sexec matrix; call loaddata; call echooff; call load(ddot :staging); y=dlog(realgdp); /; Greene(2008) page 754 c0 - c4 function ddot(n,x,ixbegin,y,iybegin); /; /; Simulate ddot functionality /; /; n => # of terms /; x => vector 1 /; ixbegin => x term to begin with /; y => vector 1 /; iybegin => y term to begin with /; /; Built 11 October 2008 /; if((ixbegin+n-1).gt.norows(x))then; call epprint('ERROR: in ddot n+ixbegin > norows(x)'); dd=missing(); go to endd; endif; if((iybegin+n-1).gt.norows(y))then; call epprint('ERROR: in ddot n+iybegin > norows(y)'); dd=missing(); go to endd; endif; if(ixbegin.lt.1.or.iybegin.lt.1)then; call epprint('ERROR: ixbegin and iybegin must be gt 1'); go to endd; endif; dd=vfam(x(integers(ixbegin,ixbegin+n-1))) * vfam(y(integers(iybegin,iybegin+n-1))) ; endd continue; return(dd); end; subroutine geteq(y); yy=y; dy=dif(yy); tt=dfloat(integers(1,norows(yy)-1)); yy(1)=missing(); yy=goodrow(yy); call olsq(yy tt yy{1} dy{1} :print); r=%res; c0= ddot(norows(r),r,1,r,1)/dfloat(norows(r)); c=array(4:); do i=1,4; c(i)=ddot(norows(r)-i,r,1,r,1+i)/dfloat(norows(r)); enddo; call print(c0,c); return; end; call geteq(y); b34srun; == ==DF_PP Dickey Fuller, Phillips Perron and KPSS Tests /; Testing DF and PP using Hamilton as a reference %b34slet runsas = 0; b34sexec options ginclude('b34sdata.mac') member (hamilton1); b34srun; b34sexec matrix; call loaddata; call echooff; call load(df_pp :staging); call load(kpss); call load(df_gls); n=15; call df_pp(tbill, 'Tbill Rate Data From Hamilton',n); call print(' ':); call df_pp(ln_gnp,'log(gnp) Data From Hamilton',n); b34srun; %b34sif(&runsas.ne.0)%then; b34sexec options open('testsas.sas') unit(29) disp=unknown$ b34srun$ b34sexec options clean(29) $ b34seend$ b34sexec pgmcall idata=29 icntrl=29$ sas $ * sas commands next ; pgmcards$ proc arima; identify var=tbill stationarity=(pp=4); run; proc arima; identify var=ln_gnp stationarity=(pp=4); run; proc arima; identify var=tbill stationarity=(dicky=4); run; proc arima; identify var=ln_gnp stationarity=(dicky=4); run; b34sreturn$ b34srun $ b34sexec options close(29)$ b34srun$ /$ the next card has to be modified to point to sas location /$ be sure and wait until sas gets done before letting b34s resume /$ *************************************************************** b34sexec options dodos('start /w /r sas testsas' ) dounix('sas testsas' ) $ b34srun$ b34sexec options npageout noheader writeout(' ','output from sas',' ',' ') writelog(' ','output from sas',' ',' ') copyfout('testsas.lst') copyflog('testsas.log') dodos('erase testsas.sas','erase testsas.lst','erase testsas.log') dounix('rm testsas.sas','rm testsas.lst','rm testsas.log') $ b34srun$ b34sexec options header$ b34srun$ %b34sendif; == ==DFVALUES Dickey-Fuller Critical Values b34sexec matrix; call load(dfvalues :staging); call echooff; /; Uses RATS logic. Values not 100% exact but close! do i=3,1,-1; n1=14; n2=100; cvalue1 =array(n2-n1+1:); cvalue5 =array(n2-n1+1:); cvalue10=array(n2-n1+1:); noob =array(n2-n1+1:); ihave=1; do j=n1,n2; call dfvalues(i,cvalue,j); cvalue1(ihave) =cvalue(1); cvalue5(ihave) =cvalue(2); cvalue10(ihave)=cvalue(3); noob(ihave)=j; ihave=ihave+1; enddo; if(i.eq.1)call print('Critical Values for Intercept and Trend':); if(i.eq.2)call print('Critical Values for Intercept':); if(i.eq.3)call print('Critical Values for No Intercept':); call tabulate(noob,cvalue1,cvalue5,cvalue10); enddo; b34srun; == ==DISPMARS Display Mars Output b34sexec options ginclude('b34sdata.mac') member(gas); b34srun; b34sexec matrix; call loaddata; call load(dispmars :staging); call echooff; call mars(gasout gasin{1 to 6} gasout{1 to 6} :nk 60 :mi 3 :print); call dispmars; b34srun; == ==DISPMARS2 Display Mars where cat Variables are Used b34sexec matrix; /; /; Test file to illustrate Mars / CAT Variable Processing. /; /; Internal list of contrl variables provided to help in user /; program writting /; call load(dispmars :staging); call echooff; n=30; x=rn(array(n:)); d= array(n:); where(x.gt.1.0)d=1.; d=d+1.; y=50.+x+d; /; Perfect fit continuous case 0-1 call olsq(y x d :print); yhat1=%yhat; call mars(y x d :print :nk 20 :mi 2 ); call dispmars; call tabulate(%names,%typevar,%states); yhat2=%yhat; /; Same model solved as a cat model 0-1 call character(lx 'd 0 catnores'); call mars(y x d :print :nk 20 :mi 2 :lx lx :nmcv 2 :ntcv 2); call dispmars; call tabulate(%names,%typevar,%states %coef %knot %typek,%parent,%varink); yhat3=%yhat; /; Test to get yhat3 mask1=(x .ge. %knot(1)); mask2=(x .lt. %knot(1)); mask3=(d .eq. %states(1)); term1=x-%knot(1); term2=%knot(1)-x; yhat3a=%coef(1)+%coef(2)*mask1*term1+ %coef(3)*mask2*term2 + %coef(4)*mask3; diff=yhat3-yhat3a; call tabulate(y yhat1 yhat3 yhat3a diff x d mask1 mask2 mask3); call tabulate(y x d yhat1,yhat2,yhat3); call print('#################################################':); call print('Three State Problem Test Case ':); call print('#################################################':); d= array(n:); where(x.ne.0.)d=1.; where(x.lt.0.0)dm=-50.; where(x.gt. .5)dd=20.; d=8.*d+100.*dd +77.*dm; y=50.+x+d; call olsq(y x d :print); yhat1=%yhat; call mars(y x d :print :nk 20 :mi 2 ); call dispmars; call tabulate(%names,%typevar,%states %coef %knot %typek,%parent,%varink); yhat2=%yhat; call character(lx 'd 0 catnores'); call mars(y x d :print :nk 20 :mi 2 :lx lx :nmcv 3 :ntcv 3); call dispmars; call tabulate(%names,%typevar,%states %coef,%knot, %typek,%parent,%varink); call print(%cknot); /; Test to get yhat3 yhat3=%yhat; /$ call names; mask1 =(d .eq. %states(1)); mask2 =(d .eq. %states(3)); mask3 =(x .gt. %knot(5) ); mask4 =(x .lt. %knot(5) ); term2 =x-%knot(5); term3 =%knot(6)-x; yhat3a=%coef(1)+%coef(2)*mask1 +%coef(4)*mask2 +%coef(6)*mask3*term2 +%coef(7)*mask4*term3; diff=yhat3-yhat3a; call tabulate(mask1,mask2,mask3,mask4,d,x,yhat3,yhat3a,diff); yhat3=%yhat; b34srun; == ==DISPMARS3 Power Problem b34sexec options ginclude('b34sdata.mac') member(power); b34srun; b34sexec matrix; call loaddata; call load(dispmars :staging); call print(dispmars); call print(disp_max); call print(getcat); call echooff; call print('*******************************************************'); call print('** Analysis Performed on Variable: DAYLOAD'); call print('*******************************************************'); /$ **************************************************************$/ /$ Set span for sample: DAYLOAD /$ **************************************************************$/ iorigins=integers(1, 730); DAYLOAD = DAYLOAD(iorigins); DATE = DATE(iorigins); DOW = DOW(iorigins); RECRDS = RECRDS(iorigins); COOLDGRE = COOLDGRE(iorigins); HEATDGRE = HEATDGRE(iorigins); TEMPERTR = TEMPERTR(iorigins); /$ ******************************************************* /$ _holdout is the number of forecasts /$ _maxlag is the longest lag in model /$ _nlags is number of lags for diagnostics /$ _knots is the maximum number of basis function /$ _mi is the maximum number of interactions /$ _ms is the minimum span between each knot /$ _ffmars holds the MARS forecasts /$ _ffols holds the OLS forecasts for comparison /$ _actuals holds the actuals for comparison /$ ******************************************************* _holdout=90; _nlags=12; _maxlag=1; _knots=120; _mi=2; _ms=0; nn=vector(_holdout+1:); _ffmars=vector(_holdout+1:); _ffols=vector(_holdout+1:); /$ ******************************************************** /$ dload1 is a copy of the original series /$ temp1 is a copy of the original series /$ ******************************************************** dload1=dayload; temp1=tempertr; /$ ******************************************************** /$ Create Day of Week Dummies /$ ******************************************************** DTMP=vector(7: 1 1 1 1 1 1 1); DMAT=diagmat(DTMP); DMAT1=DMAT; n=norows(dload1)/7; do i=1,n; DMAT=catrow(DMAT,DMAT1); enddo; D1=DMAT(,1); D2=DMAT(,2); D3=DMAT(,3); D4=DMAT(,4); D5=DMAT(,5); D6=DMAT(,6); d6(2)=2.; itmp=integers(norows(DAYLOAD)-_holdout+1,norows(DAYLOAD)); _actuals=vfam(DAYLOAD(itmp)); jj=-1*_holdout; icount=1; /$ ******************************************************** /$ Specify model components using character/argument /$ ******************************************************** call character(mymodel, 'TEMPERTR{0 to 1} DAYLOAD{1 } '); call character(dummies, 'D1{0} D2{0} D3{0} D4{0} D5{0} D6{0}'); call character(dummies2, 'DOW'); call character(lx, 'd1 0 catnores' 'd2 0 catnores' 'd3 0 catnores' 'd4 0 catnores' 'd5 0 catnores' 'd6 0 catnores'); call character(lx2, 'dow 0 catnores'); call olsq(DAYLOAD argument(mymodel) argument(dummies) :print :holdout _holdout); _olscoef=%coef; /$ ******************************************************* /$ Specify MARS model /$ ******************************************************* call mars(DAYLOAD argument(mymodel) argument(dummies) :print :nk _knots :mi _mi :ms _ms ); call dispmars; call print('at one *********************************':); call tabulate(%names,%typevar,%states %coef,%knot, %typek,%parent,%varink); call mars(DAYLOAD argument(mymodel) argument(dummies) :print :nk _knots :mi _mi :ms _ms :lx lx :nmcv 3 :ntcv 13 :isetfm 170000); call dispmars; call print('at two *********************************'); call tabulate(%names,%typevar,%states %coef,%knot, %typek,%parent,%varink); call mars(DAYLOAD argument(mymodel) argument(dummies2) :print :nk _knots :mi _mi :ms _ms :lx lx2 :nmcv 7 :ntcv 7 :isetfm 170000 ); call dispmars; call print('at three *******************************'); call tabulate(%names,%typevar,%states %coef,%knot, %typek,%parent,%varink); b34srun; == ==F_MARSLG Fixes MARSPLINE 0-1 output /; /; Shows how MARS logit and MARSPLINE on a catagolical left hand /; variable correlate with actual data and with each other. /; The tentative conclusion is that MARSPLINE using the adjustment will /; outperform OLS. More research is needed on this before "automatic" /; adjustments are made. A tentative version of f_marslg is in the /; staging library but is subject to major changes. The driving job /; is f_marslg. The complete job is marspline8 /; b34sexec options ginclude('b34sdata.mac') macro(murder)$ b34seend$ /; b34sexec probit; model d1 = t y lf nw; b34srun; /; b34sexec loglin; model d1 = t y lf nw; b34srun; b34sexec matrix; call loaddata; call load(f_marslg :staging); yvar=afam(d1); /; this sets model call character(cc,'t y lf nw'); call olsq(yvar argument(cc) :print); %yhatols=%yhat; call probit(yvar argument(cc):print); %yhat_pb=%yhat; call mars(yvar argument(cc) :logit :nk 40 :mi 2 :print); %yhat1=%yhat; call marspline(yvar argument(cc) :nk 20 :mi 2 :df 2. :print); %yhat2=%yhat; call f_marslg(%yhat2,%yhat3,%yhat4); call tabulate(%y,%yhatols,%yhat_pb,%yhat1,%yhat2,%yhat3,%yhat4); data=catcol( %y,%yhatols,%yhat_pb,%yhat1,%yhat2,%yhat3,%yhat4); call print( 'Correlation %y,%yhatols,%yhat_pb,%yhat1,%yhat2,%yhat3,%yhat4', ccf(data)); call graph( %yhatols,%yhat_pb,%yhat1,%yhat2,%yhat3,%yhat4:nolabel); b34srun; == ==G4TEST Pena-Slate Global Validation Test b34sexec options ginclude('b34sdata.mac') member(milage); b34srun; b34sexec matrix; call loaddata; call load(g4test :staging); call echooff; call olsq(n_gal miles_lf daysbetw :print :savex ); call g4test; b34srun; b34sexec options ginclude('b34sdata.mac') member(salinity); b34srun; b34sexec matrix; call loaddata; call load(g4test :staging); call echooff; call olsq(salin lagsalin trend h2oflow :print :savex); call g4test; b34srun; b34sexec options ginclude('b34sdata.mac') member(wine); b34srun; b34sexec matrix; call loaddata; call load(g4test :staging); call echooff; call olsq(h_a_rate consump :print :savex); call g4test; b34srun; b34sexec options ginclude('b34sdata.mac') member(gas); b34srun; b34sexec matrix; call loaddata; call load(g4test :staging); call echooff; call olsq(gasout gasin{1 to 6} gasout{1 to 6} :print :savex); call g4test; b34srun; b34sexec options ginclude('b34sdata.mac') member(res72); b34srun; b34sexec matrix; call loaddata; call load(g4test :staging); call echooff; call olsq(lnq lnl lnk :print :savex); call g4test; call olsq(lnq lnl lnk time :print :savex); call g4test; call olsq(lnq lnl lnk time lnrm1 :print :savex); call g4test; call olsq(lnq lnl lnk time lnrm2 :print :savex); call g4test; b34srun; == ==G4TEST_1 Illustrates Leverage and RR calculations b34sexec options ginclude('b34sdata.mac') member(milage); b34srun; b34sexec matrix; call loaddata; call load(g4test :staging); call echooff; /; Note: /; call g4test2(res,y,yhat,beta,s1, s2, s3, s4,g4,x, /; ps1,ps2,ps3,ps4,pg4,iprint); call olsq(n_gal miles_lf daysbetw :print :savex ); call g4test; call g4move; call g4plot; call tabulate(%g4dif, %pg4dif, %s1dif ,%ps1dif, %s2dif, %ps2dif,%s3dif, %ps3dif,%s4dif, %ps4dif); call g4rr; call g4rrplot; call tabulate(%rg4, %rpg4, %rs1, %rps1, %rs2, %rps2, %rs3, %rps3, %rs4, %rps4); b34srun; == ==G4TEST_2 Generated Data %b34slet noob=15000$ b34sexec data noob=%b34seval(&noob)$ build y1 y2 x z e1 e2$ gen e1=rn()$ gen e2=rn()$ gen x =10*rn()$ gen z =10*rn()$ gen y1 = 10 + 5*x + 5*z + 50*e1 $ gen if(x .gt. 0) y2= 10 + 5*x + 5*z + 50*e2$ gen if(x .le. 0) y2= 10 -10*x + 5*z + 50*e2$ b34srun$ b34sexec matrix; call loaddata; call load(g4test :staging); call echooff; call olsq(y1 x z :print :savex); call g4test; call olsq(y2 x z :print :savex); call g4test; b34srun; == ==G4TEST_3 Gas Data b34sexec options ginclude('gas.b34') $ b34srun$ b34sexec matrix; call loaddata; call echooff; maxlag=6; call load(g4test :staging); call olsq(gasout gasin{1 to maxlag} gasout{1 to maxlag} :print :savex); call g4test; b34srun; == ==G4TEST_4 Friedman b34sexec options ginclude('b34sdata.mac') member(friedman); b34srun; b34sexec matrix; call loaddata; call load(g4test :staging); call echooff; call olsq(y x1 x2 x3 x4 x5 :savex :print); call g4test; b34srun; == ==G4TEST_5 Brieman/ Cheveland Devlin Nitrous Oxide Data b34sexec data noob=88 heading('nitrous oxide data') nohead filef=@@$ * data discussed page 133 in breiman (1991)$ * data came from Cleveland and Devlin (1988)$ * data transformed such that y = nox**(1/3) - median(nox**(1/3))$ input y e_ratio c_ratio$ label y = 'nitrous oxides in exhaust'$ label e_ratio = 'equivalence ratio '$ label c_ratio = 'compression ratio '$ datacards$ 3.741 .907 12. 2.295 .761 12. 1.498 1.108 12. 2.881 1.016 12. .76 1.189 12. 3.12 1.001 9. .638 1.231 9. 1.17 1.123 9. 2.358 1.042 12. .606 1.215 12. 3.669 .93 12. 1. 1.152 12. .981 1.138 15. 1.192 .601 18. .926 .696 7.5 1.59 .686 12. 1.806 1.072 12. 1.962 1.074 15. 4.028 .934 15. 3.148 .808 9. 1.836 1.071 9. 2.845 1.009 7.5 1.013 1.142 7.5 .414 1.229 18. .812 1.175 18. .374 .568 15. 3.623 .977 15. 1.869 .767 7.5 2.836 1.006 7.5 3.567 .893 9. .866 1.152 15. 1.369 .693 15. .542 1.232 15. 2.739 1.036 15. 1.2 1.125 15. 1.719 1.081 9. 3.423 .868 9. 1.634 .762 7.5 1.021 1.144 7.5 2.157 1.045 7.5 3.361 .797 18. 1.39 1.115 18. 1.947 1.07 18. .962 1.219 18. .571 .637 9. 2.219 .733 9. 1.419 .715 9. 3.519 .872 9. 1.732 .765 7.5 3.206 .878 7.5 2.471 .811 7.5 1.777 .676 15. 2.571 1.045 18. 3.952 .968 18. 3.931 .846 15. 1.587 .684 15. 1.397 .729 7.5 3.536 .911 7.5 2.202 .808 7.5 .756 1.168 7.5 1.62 .749 7.5 3.656 .892 7.5 2.964 1.002 7.5 3.76 .812 18. .672 1.23 18. 3.677 .804 18. 3.517 .813 12. 3.29 1.002 12. 1.139 .696 9. .727 1.199 9. 2.581 1.03 9. .923 .602 15. 1.527 .694 15. 3.388 .816 15. 2.085 1.037 15. .966 1.181 15. 3.488 .899 7.5 .754 1.227 7.5 .797 1.18 9. 2.064 .795 7.5 3.732 .99 18. .586 1.201 18. .561 .629 7.5 .563 .608 9. .678 .584 12. .37 .562 15. .53 .535 18. 1.9 .655 18. b34sreturn$ b34seend$ b34sexec matrix; call loaddata; call load(g4test :staging); call echooff; call olsq( y e_ratio c_ratio :print :savex); call g4test; b34srun; == ==G4TEST_6 Ozone data /; /; OLS/GAM/MARS/ACE /; b34sexec options ginclude('b34sdata.mac') member(b_ozone); b34srun; b34sexec matrix; call loaddata; call echooff; call load(g4test :staging); call olsq(ozone vh wind humidity temp ibh dpg ibt vis doy :print :savex); call g4test; b34srun; == ==G4TEST_7 Tree Data b34sexec options ginclude('b34sdata.mac') member(trees); b34srun; b34sexec matrix; call loaddata; call load(g4test :staging); call echooff; call olsq(volume girth height :print :savex); call g4test; b34srun; == ==G4TEST_8 SO4 and Latitude and Longitude /; Study of SO4 and Latitude and Longitude /; Data Studied byDong Xiang at SAS Institute /; 'Fitting Generalized Additive Models with GAM Procedure' /; SAS paper 256-26 b34sexec options ginclude('b34sdata.mac') member(gam_6); b34srun; b34sexec matrix; call loaddata; call load(g4test :staging); call echooff; call olsq(so4 latt long :print :savex); call g4test; b34srun; == ==GARCHF Forecasting a ARCH / GARCH Model b34sexec scaio readsca /; file('/usr/local/lib/b34slm/findat01.mad') file('c:\b34slm\examples\findat01.mad') dataset(m_intc); b34srun; %b34slet dob34s = 1; %b34slet dorats = 0; /$ The following Script was built by Workbench %b34sif(&dob34s.eq.1)%then; b34sexec matrix; call loaddata; call load(garchf :staging); call echooff; call print('*******************************************************'); call print('** Analysis Performed on Variable: INTEL'); call print('Options: ARCH , Diags: ON, Graphs: ON'); call print('*******************************************************'); /$ **************************************************************$/ /$ Set span for time series: INTEL /$ **************************************************************$/ iorigins=integers(1, 300); INTEL = INTEL(iorigins); /$ **************************************************************$/ /$ Use mean/variance for CONSTANT starting values in model /$ **************************************************************$/ YMean=Mean(INTEL); YSigma2=Variance(INTEL-YMean); call garchest(res1, res2, INTEL,func,3,nbad :cparms array(:YMean, YSigma2) :nma idint(array(:1 2 )) :gmaorder idint(array(:1 2 3)) :print); call print('Residual Sum of Squares is:',sumsq(goodrow(res1)):); /; Forecasting nfore=10; iprintf=1; mtype=1; nosqrt=0; call garchf(mtype,res1,res2,intel,nfore,fore1,fore2,iobs,%nparm, %cname,%coef,%corder,%nar,%nma,%ngar,%ngma,%nmu,%con, nosqrt,iprintf); /; --------------------------------------------------------------;/ /$ **************************************************************$/ /$ Set number of lags for ACF/PACF/Q-Stat /$ **************************************************************$/ _nlags=24; /$ **************************************************************$/ /$ Produce ACF and PACF statistics for INTEL /$ **************************************************************$/ call print('*******************************************************'); call print('** Diagnostics/Summary Stats for Dependent Var: INTEL'); call print('*******************************************************'); call print('Print--> Describe Dependent Var: INTEL '); call describe(INTEL :print); call print('Print--> ACF,Std.Err,PACF,Q-Stat,Prob.Q for INTEL'); acfy=acf(INTEL,_nlags,sey,pacfy,mqy,pmqy); pmqy=1.-pmqy; call tabulate(acfy,sey,pacfy,mqy,pmqy); /$ **************************************************************$/ /$ Standardize residuals of 1st moment model /$ **************************************************************$/ resa=goodrow(goodrow(res1)/dsqrt(goodrow(res2))); /$ **************************************************************$/ /$ Build LaGrange Multiplier Test for 1st moment residuals /$ **************************************************************$/ n=2; lmvalue=array(n:); lag=idint(array(n:)); prob=array(n:); do i=1,n; lag(i)=i; call lm(goodrow(resa), value,i,pp); lmvalue(i)=value; prob(i)=pp; enddo; call print('*******************************************************'); call print('** Diagnostics of standardized residuals'); call print('*******************************************************'); call print('Print--> Describe standardized resid. series'); call describe(resa :print); call print('Print--> Engle LM Test for standardized res1 series'); call tabulate(lag,lmvalue,prob); /$ **************************************************************$/ /$ Compute ACF, PACF, Q-stats for residuals in 1st moment model /$ **************************************************************$/ call print('Print--> ACF, Std.Err, PACF, Q-Stat, Prob.Q'); acfa=acf(resa,_nlags,sea,pacfa,mqa,pmqa); pmqa=1.-pmqa; call tabulate(acfa,sea,pacfa,mqa,pmqa); /$ **************************************************************$/ /$ Square standardized residuals of 1st moment model /$ **************************************************************$/ resb=resa**2.; call print('*******************************************************'); call print('** Diagnostic testing of squared standardized residuals'); call print('*******************************************************'); /$ **************************************************************$/ /$ Compute ACF, PACF, Q-stats for squared standardized residuals /$ **************************************************************$/ call print('Print--> ACF, Std.Err, PACF, Q-Stat, Prob.Q'); acfb=acf(resb,_nlags,seb,pacfb,mqb,pmqb); pmqb=1.-pmqb; call tabulate(acfb,seb,pacfb,mqb,pmqb); call print('*******************************************************'); call print('** Display graphics for Dependent Variable'); call print('*******************************************************'); call graph(INTEL :file 'yvar.wmf' :pspaceon :pgyscaleright 'i' :pgborder :pgxscaletop 'i' :nocontact :colors black bblue :heading 'Plot of INTEL'); call graph(acfy, sey :file 'acfy.wmf' :pspaceon :pgyscaleright 'i' :pgborder :pgxscaletop 'i' :histscale integers(0,_nlags,2) :overlay acfplot :colors black bblue bred :heading 'ACF of INTEL'); call graph(pacfy, sey :file 'pacfy.wmf' :pspaceon :pgyscaleright 'i' :pgborder :pgxscaletop 'i' :histscale integers(0,_nlags,2) :overlay acfplot :colors black bblue bred :heading 'Partial ACF of INTEL'); call print('*******************************************************'); call print('** Display graphics for 1st Moment Residuals'); call print('*******************************************************'); call graph(resa :file 'resa.wmf' :pspaceon :pgyscaleright 'i' :pgborder :pgxscaletop 'i' :nocontact :colors black bblue :heading 'Plot of 1st Moment Residuals'); call graph(acfa, sea :file 'acfa.wmf' :pspaceon :pgyscaleright 'i' :pgborder :pgxscaletop 'i' :histscale integers(0,_nlags,2) :overlay acfplot :colors black bblue bred :heading 'ACF of 1st Moment Residuals'); call graph(pacfa, sea :file 'pacfa.wmf' :pspaceon :pgyscaleright 'i' :pgborder :pgxscaletop 'i' :histscale integers(0,_nlags,2) :overlay acfplot :colors black bblue bred :heading 'PACF of 1st Moment Residuals'); call graph(mqa :file 'mqa.wmf' :pspaceon :pgyscaleright 'i' :pgborder :pgxscaletop 'i' :nocontact :colors black bblue :heading 'Q-Stats from 1st Moment Residuals'); call print('*******************************************************'); call print('** Display graphs for Squared Standardized Residuals'); call print('*******************************************************'); call graph(resb :file 'resb.wmf' :pspaceon :pgyscaleright 'i' :pgborder :pgxscaletop 'i' :nocontact :colors black bblue :heading 'Plot of Squared Standardized Residuals'); call graph(acfb, seb :file 'acfb.wmf' :pspaceon :pgyscaleright 'i' :pgborder :pgxscaletop 'i' :histscale integers(0,_nlags,2) :overlay acfplot :colors black bblue bred :heading 'ACF of Squared Standardized Residuals'); call graph(pacfb, seb :file 'pacfb.wmf' :pspaceon :pgyscaleright 'i' :pgborder :pgxscaletop 'i' :histscale integers(0,_nlags,2) :overlay acfplot :colors black bblue bred :heading 'Partial ACF of Squared Standardized Residuals'); call graph(mqb :file 'mqb.wmf' :pspaceon :pgyscaleright 'i' :pgborder :pgxscaletop 'i' :nocontact :colors black bblue :heading 'Q-Stats from Squared Standardized Residuals'); b34srun; %b34sendif; %b34sif(&dorats.eq.1)%then; b34sexec options open('rats.dat') unit(28) disp=unknown$ b34srun$ b34sexec options open('rats.in') unit(29) disp=unknown$ b34srun$ b34sexec options clean(28)$ b34srun$ b34sexec options clean(29)$ b34srun$ b34sexec pgmcall$ rats passasts PCOMMENTS('* ', '* Data passed from B34S(r) system to RATS', '* ') $ PGMCARDS$ ************************ * This program performs ARCH(3) Estimation. ************************ * all 0 300:1 * open data m-intc.dat * * data(org=obs) / rt dat * * define the log-likelihood recursively * set rt = intel set h = 0.0 * nonlin mu a0 a1 a2 a3 * frml at = rt(t) - mu frml gvar = a0 + a1*at(t-1)**2 + a2*at(t-2)**2 + a3*at(t-3)**2 frml garchlog = -0.5*log(h(t)=gvar(t))-0.5*at(t)**2/h(t) * smpl 4 300 * * initial values * compute mu = 0.01, a0 = 0.01, a1 = 0.2, a2 = 0.1, a3 = 0.1 * maximize(method=bhhh,recursive,iterations=100) garchlog * * Select the estimation method: bhhh * Specify that the function is recursively defined and increase the number * of iterations from 20 (default) to 100. * set et = at(t) set fvar = gvar(t) set resi = at(t)/sqrt(fvar(t)) set resisq = resi(t)*resi(t) *** Checking standardized residuals *** cor(qstats,number=20,span=10) resi *** Checking squared standardized residuals *** cor(qstats,number=20,span=10) resisq *** Last few observations needed for forecasts *** print 295 300 rt et resi fvar b34sreturn$ b34srun $ b34sexec options close(28)$ b34srun$ b34sexec options close(29)$ b34srun$ b34sexec options /$ dodos(' rats386 rats.in rats.out ') dodos('start /w /r rats32s rats.in /run') dounix('rats rats.in rats.out')$ B34SRUN$ b34sexec options npageout WRITEOUT('Output from RATS',' ',' ') COPYFOUT('rats.out') dodos('ERASE rats.in','ERASE rats.out','ERASE rats.dat') dounix('rm rats.in','rm rats.out','rm rats.dat') $ B34SRUN$ %b34sendif; == ==GARCH2PA Test Problem auto 2 pass Garch b34sexec options ginclude('gas.b34'); b34srun; /; /; See also AUTOGARCH template /; b34sexec matrix; call loaddata; call echooff; call load(garch2pa :staging); /; /; subroutine arguments and help notes /; /; subroutine garch2pa(x,seas,nac,noprint,nf,nt,res1,ressq1,res2); /; /; Automatic 2 Pass GARCH Modeling /; /; Note: Most options not enabled. See autogarch in matrix menus. /; /; x - Series input /; seas - Seasonal. Must be integer ge 0 input /; nac - # Autocorrelations input /; noprint - =0 for no print input /; nf - # of forecasts. Default = 1 input /; nt - Forecast time base. =0 => last ops input /; Note: second moment foirecast done at period end /; If no print in effect, no forecasts done /; res1 - First Monent Residual output /; ressq1 - First Moment Residual squared output /; res2 - Second Moment residual output /; /; Routine built 4 July 2004 /; ***************************************************************** iseas=0; nac=24; nprint=1; nf=10; nt=0; call garch2pa(gasout,iseas,nac,nprint,nf,nt,res1,ressq1,res2); b34srun; == ==GENGARCH Simulate GARCH & ARCH Models B34SEXEC MATRIX DISPLAY=COL80MEDIUM; /$ /$ Define the number of observations to be simulated /$ in the parameter N. /$ N=1010; call load(gengarch :staging); s2 =5.0; idrop=10; gar =array(:.6 ); gma =array(:.3); call echooff; call gengarch(n,s2,gar,gma,idrop,e,h); /$ **************************************************************$/ /$ Specify and estimate GARCH model to test simulation /$ **************************************************************$/ call graph(e :heading 'Vector of GARCH(p,q) innovations'); call graph(h :heading 'Associated conditional variance process'); c1=5.; c2=.5; res1=array(norows(e):); res2=array(norows(e):); call garchest(res1, res2, e,func,2,nbad :garorder idint(array(:1)) :gmaorder idint(array(:1)) :cparms array(:c1,c2) :print); B34SRUN; == ==GENGARCH2 GARCH(2,1) Simulation B34SEXEC MATRIX DISPLAY=COL80MEDIUM; /$ /$ Define the number of observations to be simulated /$ in the parameter N. /$ /$ N=10000; N=1000; call load(gengarch :staging); s2 =5.0; idrop=10; gar =array(:.4,0.1); /$gar =array(: .6); gma =array(:0.2); call echooff; call gengarch(n,s2,gar,gma,idrop,e,h); /$ **************************************************************$/ /$ Specify and estimate GARCH model to test simulation /$ **************************************************************$/ call graph(e :heading 'Vector of GARCH(p,q) innovations'); call graph(h :heading 'Associated conditional variance process'); res1=array(norows(e):); res2=array(norows(e):); call garchest(res1, res2, e,func,2,nbad :garorder idint(array(:1 2)) :gmaorder idint(array(:1)) :cparms array(:c1,c2) :print); b34srun; == ==GETCHAR Unpack a string into a Character Array b34sexec matrix; call load(ntokin :staging); call load(getr16 :staging); call load(getchar :staging); call echooff; /$ Loading Answers into Character then reading in real*16 /$ and character call character(cc,' -1467.48961422980 298.084530995537 -2772.17959193342 559.779865474950 -2316.37108160893 466.477572127796 -1127.97394098372 227.204274477751 -354.478233703349 71.6478660875927 -75.1242017393757 15.2897178747400 -10.8753180355343 2.23691159816033 -1.06221498588947 0.221624321934227 -0.670191154593408E-01 0.142363763154724E-01 -0.246781078275479E-02 0.535617408889821E-03 -0.402962525080404E-04 0.896632837373868E-05' ); call ntokin(cc,nfind,0,ibad); call print(nfind ); call getr16(cc,nfind,ans,ibad); call getchar(cc,nfind,cx,IBAD); do i=1,norows(ans); call fprint(:clear :col 4 :string 'Real*16 ' :col 20 :display ans(i) '(e25.16)' :print :clear :col 4 :string 'Character Array' :col 23 :string cx(i,) :print); enddo; b34srun; == ==GETDATA Simple real*8 data read example b34sexec matrix; * Simple read where a matrix on 10 by 200 is built and read back; * Data is real*8; * Job shows how to beat 98 variable limit; * Note: Before reading, structure of object must be known!!!! ; call echooff; nobs=10; nvar=200; test=rn(array(nobs,nvar:)); call open(70,'testdata'); call rewind(70); call write(test,70); call print(test); /; Reading data call load(getdata: staging); file='testdata'; nrows=nobs; ncols=nvar; call getdata(file,nrows,ncols,x); call print(x); b34srun; == ==GETDATA2 Simple real*16 data read example b34sexec matrix; * Data is real*16; * Note: Before reading, structure of object must be known!!!! ; call echooff; nobs=10; nvar=8; test=r8tor16(rn(array(nobs,nvar:))); call open(70,'testdata'); call rewind(70); call write(test,70); call print(test); /; Reading data call load(getdata2: staging); file='testdata'; nrows=nobs; ncols=nvar; call getdata(file,nrows,ncols,x); call print(x); b34srun; == ==GETDATA2B Compairs two ways to read real*16 Data b34sexec matrix; * Data is real*16; * Note: Before reading, structure of object must be known!!!! ; call echooff; nobs=10; nvar=8; test=r8tor16(rn(array(nobs,nvar:))); /; test=rn(array(nobs,nvar:)); call open(70,'testdata'); call rewind(70); call write(test,70); call print(test); /; Reading data call load(getdata2: staging); file='testdata'; nrows=nobs; ncols=nvar; call getdata2(file,nrows,ncols,x); call print(x); /$ /$ Reads a character array into real*16 and real*8. Tests input /$ call load(ntokin :staging); call load(getr16 :staging); call load(getr8 :staging); /; call echooff; cc=c1array(5000:); call open(70,'testdata'); call rewind(70); call read(cc,70,'(5000a1)'); call print(cc); call ntokin(cc,nfind,0,ibad); call print('All done tokenizing'); call print('# found was ',nfind); call getr16(cc,nfind,x16,ibad); call print(cc); call print(x16); call getr8(cc,nfind,x8,ibad); call print(x8); /; Look at Difference call print('Difference between two reads'); /; Repack note how loaded x16=array(nocols(x),norows(x):x16); call print((transpose(x)-x16)); b34srun; == ==GETR16 Unpack a string into Real*16 Array /$ /$ Reads a character array into real*16 and real*8. Tests input /$ b34sexec matrix; call character(cc,' 1 0 6 1 63 2 364 3 1365 4 3906 5 9331 6 19608 7 37449 8 66430 9 111111 10 177156 11 271453 12 402234 13 579195 14 813616 15 1118481 16 1508598 17 2000719 18 2613660 19 3368421 20'); call print(cc); call load(ntokin :staging); call load(getr16 :staging); call load(getr8 :staging); call echooff; call ntokin(cc,nfind,0,ibad); call print('All done tokenizing'); call print('# found was ',nfind); call getr16(cc,nfind+10,x16,ibad); call print(cc); call print(x16); call getr8(cc,nfind,x8,ibad); call print(x8); /$ repack xm=matrix(21,2:goodrow(x16)); call print(xm); b34srun; == ==GETR8 Unpack a string into Real*8 Array /$ /$ Reads a character array into real*16 and real*8. Tests input /$ b34sexec matrix; call character(cc,' 1 0 6 1 63 2 364 3 1365 4 3906 5 9331 6 19608 7 37449 8 66430 9 111111 10 177156 11 271453 12 402234 13 579195 14 813616 15 1118481 16 1508598 17 2000719 18 2613660 19 3368421 20'); call print(cc); call load(ntokin :staging); call load(getr16 :staging); call load(getr8 :staging); call echooff; call ntokin(cc,nfind,0,ibad); call print('All done tokenizing'); call print('# found was ',nfind); call getr16(cc,nfind+10,x16,ibad); call print(cc); call print(x16); call getr8(cc,nfind,x8,ibad); call print(x8); /$ repack xm=matrix(21,2:goodrow(x16)); call print(xm); b34srun; == ==GETVPA Unpack a string into VPA Array /$ /$ Reads a character array into VPA. Tests input /$ b34sexec matrix; call character(cc,' 1 0 6 1 63 2 364 3 1365 4 3906 5 9331 6 19608 7 37449 8 66430 9 111111 10 177156 11 271453 12 402234 13 579195 14 813616 15 1118481 16 1508598 17 2000719 18 2613660 19 3368421 20'); call print(cc); call load(ntokin :staging); call load(getvpa :staging); call load(getr16 :staging); call echooff; call ntokin(cc,nfind,0,ibad); call print('All done tokenizing'); call print('# found was ',nfind); call getr16(cc,nfind,x16,ibad); call print(cc); call print(x16); call getvpa(cc,nfind,vpa_x,ibad); call print(vpa_x); /$ repack xm=matrix(nfind/2,2:x16); call print(xm); vpa_xm=matrix(nfind/2,2:vpa_x); call print(vpa_xm); b34srun; == ==GETVPA Very Hard OLS Data. /; /; Problem very very hard!!! /; Illustrates reading into real*8 and then into real*16 and VPA /; Reading from character into real*16 /; Reading from chapter into VPA /; b34sexec data heading('Bruce McCullough Test Case'); input y x; datacards; 10000000001 1000000000.000 10000000002 1000000000.000 10000000003 1000000000.900 10000000004 1000000001.100 10000000005 1000000001.010 10000000006 1000000000.990 10000000007 1000000001.100 10000000008 1000000000.999 10000000009 1000000000.000 10000000010 1000000000.001 b34sreturn; b34srun; b34sexec matrix; call loaddata; call echooff; /; /; Data reread into character to see what happens /; datacards; 10000000001 1000000000.000 10000000002 1000000000.000 10000000003 1000000000.900 10000000004 1000000001.100 10000000005 1000000001.010 10000000006 1000000000.990 10000000007 1000000001.100 10000000008 1000000000.999 10000000009 1000000000.000 10000000010 1000000000.001 b34sreturn; y16=r8tor16(y); x16=r8tor16(x); /; Cholesky on real*8 fails /; call olsq(y x :print); /; QR fails /; call olsq(y x :qr :print); /; real*16 Cholesky fails /; call olsq(y16 x16 :print); call olsq(y16 x16 :qr :print :savex); coef8_16=rollright(%coef); /; Go to vpa 64 digits and 20 digits to see if makes a difference y_vpa=vpa(%y); x_vpa=vpa(%x); beta_vpa=inv(transpose(x_vpa)*x_vpa)*transpose(x_vpa)*y_vpa; call print('Beta from VPA - default 64 digits accuracy ',beta_vpa); call print(beta_vpa(1)); call print(beta_vpa(2)); call vpaset(:ndigits 200); y_vpa=vpa(%y); x_vpa=vpa(%x); beta_vpa=inv(transpose(x_vpa)*x_vpa)*transpose(x_vpa)*y_vpa; call print('Beta from VPA - 200 digits accuracy ',beta_vpa); call print(beta_vpa(1)); call print(beta_vpa(2)); /; Reading data as character into VPA and real*16 /; Makes a difference !!!!!!!!!!!!!!!!!!!!!!!!!!! data16=r8tor16(array(20:)); call read(data16,4); call load(ntokin :staging); call load(getr16 :staging); call load(getvpa :staging); call load(getchar :staging); call vpaset(:ndigits 200); vpdata=vpa(data16); vpx=vpdata*vpa('0.0'); call rewind(4); cc=c1array(10000:); call read(cc,4); call ntokin(cc,nfind,0,ibad); call print(nfind ); call getvpa(cc,nfind,vpx,ibad); call free(cc); xm=matrix(10,2:data16); y16=array(:xm(,1)); x16=matrix(10,2:data16); x16(,1)=real16('1.'); /; This will not solve due to condition /; call olsq(y16 x16 :print :noint); /; call print('Here data read into real*16 directly - Note difference':); call olsq(y16 x16 :qr :print :noint); call print(' ':); call print('Display More digits to see what real*16 has':); call print(' ':); call fprint(:clear :display %coef(1) '(1pe48.32)' :print); call fprint(:clear :display %coef(2) '(1pe48.32)' :print); call print(' ':); call print('Compare with reading into real*8 then to real*16':); call print(' ':); call fprint(:clear :display coef8_16(1) '(1pe48.32)' :print); call fprint(:clear :display coef8_16(2) '(1pe48.32)' :print); call print(' ':); /; now run with vpa vpxm=matrix(10,2:vpx); vpy= vector(:vpxm(,1)); vpx=vpxm; vpx(,1)=vpa('1.0'); betavpa=inv(transpose(vpx)*vpx)*transpose(vpx)*vpy; call print(' ':); call print('Beta from VPA where character read into VPA':); call print('This is the right way to attack the problem':); call print('For gain compare with real*16 where character read':); call print(betavpa(1)); call print(betavpa(2)); call print(' ':); call print('Beta from VPA where vpa converted from real*8':); beta_vpa=rollright(beta_vpa);; call print(beta_vpa(1)); call print(beta_vpa(2)); call print(' ':); call print('Difference of coef # 1 and # 2':); call print('Shows effect of not reading in correct precision':); diffcoef=vpa(vector(2:)); diffcoef(1)=betavpa(1)-beta_vpa(1); diffcoef(2)=betavpa(2)-beta_vpa(2); call print(diffcoef(1)); call print(diffcoef(2)); b34srun; == ==GLM_INFO Selecting best GLM model /; /; Ridge regression also shown /; /; For further detail see Stokes (200x, Chapter 10) /; /; %b34slet dorats =0; %b34slet doridge=0; %b34slet dolasso=0; b34sexec options ginclude('gas.b34'); b34srun; b34sexec matrix; call loaddata; call load(ridge :staging); call load(lasso :staging); call load(glm_info :staging); call echooff; /; logic works for holdout > 0 k=10; holdout=100; maxlag=k; /; at issue is how are yhat values calculated within the sample treated? /; Fore purposse of this analysis they are removed call olsq(gasout gasout{1 to k} gasin{1 to k} :print :savex :holdout holdout); %yfuture=gasout(integers(norows(%x)+maxlag,norows(gasout))); olsf=mfam(%xfuture)*vfam(%coef); olsfss=sumsq(afam(%yfuture)-afam(olsf)); call print('ols forecast error sum sq ',olsfss:); nlam=100; lam_min=.000003; thr =.1e-6 ; parm = .5; ne = 20; call glm(gasout gasout{1 to k} gasin{1 to k} :print :savex :lamdamin lam_min :nlam nlam :holdout holdout :thr thr :parm parm :ne ne); call print(' ':); %yfuture=gasout(integers(norows(%x)+maxlag,norows(gasout))); call glm_info(%yfuture,%xfuture,%coef,%a0,%alm,glmf,fss,mod,1); call print('glm forecast sum of squares ',fss:); res_ols=vfam(%yfuture)-olsf; res_glm=vfam(%yfuture)-glmf; call graph( %yfuture,olsf,glmf :nolabel :nocontact :pgborder :heading 'OLS vs GLM out of sample'); call graph( res_ols res_glm :nolabel :nocontact :pgborder :heading 'OLS vs GLM residual out of sample'); call tabulate(%yfuture,olsf,glmf res_ols res_glm); /; tests on reduction loss rss_glm=array(2*k:); ne_used=array(2*k:); icount=1; do ne=1,2*k; ii=ne; call glm(gasout gasout{1 to k} gasin{1 to k} :savex :lamdamin lam_min :nlam nlam :holdout holdout :thr thr :parm parm :ne ii :print); call print(' ':); %yfuture=gasout(integers(norows(%x)+maxlag,norows(gasout))); call glm_info(%yfuture,%xfuture,%coef,%a0,%alm,glmf,fss,mod,0); rss_glm(icount)=sfam(fss); ne_used(icount)=sfam(dfloat(ii)); icount=icount+1; enddo; rss_glm =dropfirst(rss_glm,1); ne_used =dropfirst(ne_used,1); call tabulate(rss_glm,ne_used); call graph(ne_used rss_glm :grid :heading 'Loss in RSS due to model reduction' :plottype xyplot :nocontact :nolabel :pgborder); %b34sif(&doridge.ne.0)%then; lamda=2.; lamda=%alm(%lmu); call print('Lamda assumed to be ',lamda); call ridge(%y,%x,lamda,%coef,%names,%lag,ridge_c,0); call ridge(%y,%x,lamda,%coef,%names,%lag,ridge_c,1); %b34sendif; %b34sif(&dolasso.ne.0)%then; /; Redefine %coef call olsq(gasout gasout{1 to k} gasin{1 to k} :savex :holdout holdout); lamda=10.*sum(%coef)/2.0; call echooff; call lasso(%y,%x,%coef,%lcoef1,%l_t1,lamda,lresid1,1); call lasso(%y,%x,%coef,%lcoef2,%l_t2,lamda,lresid2,3); call tabulate(%names,%lag,%coef,%se,%t,%lcoef1,%l_t1,%lcoef2,%l_t2); %b34sendif; b34srun; %b34sif(&dorats.ne.0)then; b34sexec options open('rats.dat') unit(28) disp=unknown$ b34srun$ b34sexec options open('rats.in') unit(29) disp=unknown$ b34srun$ b34sexec options clean(28)$ b34srun$ b34sexec options clean(29)$ b34srun$ /; Uses logic from RATS User's Guide Version 6.1 Page 192 b34sexec pgmcall$ rats passasts pcomments('* ', '* Data passed from B34S(r) system to RATS', '* ', "display @1 %dateandtime() @33 ' Rats Version ' %ratsversion()" '* ') $ PGMCARDS$ * * Non centered results * cmoment # constant gasout{1 to 6} gasin{1 to 6} gasout linreg(cmoment) gasout # constant gasout{1 to 6} gasin{1 to 6} do row=1,13 compute %cmom(row,row)=%cmom(row,row)+2 end do row linreg(cmoment) gasout # constant gasout{1 to 6} gasin{1 to 6} b34sreturn$ b34srun $ b34sexec options close(28)$ b34srun$ b34sexec options close(29)$ b34srun$ b34sexec options /$ dodos(' rats386 rats.in rats.out ') dodos('start /w /r rats32s rats.in /run') dounix('rats rats.in rats.out')$ B34SRUN$ b34sexec options npageout WRITEOUT('Output from RATS',' ',' ') COPYFOUT('rats.out') dodos('ERASE rats.in','ERASE rats.out','ERASE rats.dat') dounix('rm rats.in','rm rats.out','rm rats.dat') $ B34SRUN$ %b34sendif; == ==GLM_INFO2 Out of sample forecasting model selection /; For further detail see Stokes (200x, Chapter 10) /; b34sexec options ginclude('gas.b34'); b34srun; b34sexec matrix; call loaddata; call load(glm_info :staging); call echooff; /; logic works for holdout > 0 /; specify model lhs = 'gasout'; rhs = ' gasout{1 to maxlag} gasin{1 to maxlag}'; maxlag=20; holdout=100; nlam=30; lam_min=.000003; thr =.1e-6 ; parm = 1.; ne = 30; /; /; test reduction effect on out of sample sum squares /; res_test=1; iprint=1; iprintv=0; /; at issue is how are yhat values calculated within the sample treated? /; Fore purposse of this analysis they are removed program gml_test; if(iprint.ne.0) call olsq(argument(lhs) argument(rhs) :print :savex :holdout holdout); if(iprint.eq.0) call olsq(argument(lhs) argument(rhs) :savex :holdout holdout); temp=argument(lhs); %yfuture=temp(integers(norows(%x)+maxlag,norows(argument(lhs)))); olsf=mfam(%xfuture)*vfam(%coef); olsfss=sumsq(afam(%yfuture)-afam(olsf)); if(iprint.ne.0)call print(' ':); call print('ols out-of-sample forecast error sum squares ',olsfss:); if(iprint.ne.0) call glm(argument(lhs) argument(rhs) :print :savex :lamdamin lam_min :nlam nlam :holdout holdout :thr thr :parm parm :ne ne ); if(iprint.eq.0) call glm(argument(lhs) argument(rhs) :savex :lamdamin lam_min :nlam nlam :holdout holdout :thr thr :parm parm :ne ne ); if(iprint.ne.0)call print(' ':); temp=argument(lhs); %yfuture=temp(integers(norows(%x)+maxlag,norows(argument(lhs)))); call glm_info(%yfuture,%xfuture,%coef,%a0,%alm,glmf,fss,mod,1); best_lam =%alm(mod); call print(' ':); call print('Model Family (0.0=ridge, 1.0=lasso) 'parm:); call print('Number of Models Considered ',nlam:); call print('Model selected based on out of sample testing ',mod:); call print('glm out-of-sample forecast error sum of squares ',fss:); call print('Lamda value for Selected Model ', best_lam:); test=%coef(,mod); where(test.eq.0.0)test=missing(); test=goodrow(test); nn=norows(test); nn2=norows(%coef(,mod)); call print('Number of Nonzero Coefficients ',nn :); call print('Number of Coefficients ',nn2:); call print(%coef(,mod)); res_ols=vfam(%yfuture)-olsf; res_glm=vfam(%yfuture)-glmf; call graph( %yfuture,olsf,glmf :nolabel :nocontact :pgborder :heading 'OLS vs GLM out of sample'); call graph( res_ols res_glm :nolabel :nocontact :pgborder :heading 'OLS vs GLM residual out of sample'); call tabulate(%yfuture,olsf,glmf res_ols res_glm); /; tests on reduction loss if(res_test.ne.0)then; rss_mod =array(ne:); rss_glm =array(ne:); ne_used =array(ne:); best_mod =array(ne:); icount=1; do ii=1,ne; if(iprintv.eq.0) call glm(argument(lhs) argument(rhs) :savex :lamdamin lam_min :nlam nlam :holdout holdout :thr thr :parm parm :ne ii ); if(iprintv.ne.0)then; call glm(argument(lhs) argument(rhs) :savex :lamdamin lam_min :nlam nlam :holdout holdout :thr thr :parm parm :ne ii :print); call print(' ':); endif; temp=argument(lhs); %yfuture=temp(integers(norows(%x)+maxlag,norows(argument(lhs)))); call glm_info(%yfuture,%xfuture,%coef,%a0,%alm,glmf,fss,mod,0); rss_glm(icount)=sfam(fss); ne_used(icount)=sfam(dfloat(ii)); best_mod(icount)=mod; rss_mod(icount)=%rss(mod); icount=icount+1; enddo; rss_glm =dropfirst(rss_glm, 1); rss_mod =dropfirst(rss_mod, 1); ne_used =dropfirst(ne_used, 1); best_mod =idint(dropfirst(dfloat(best_mod),1)); if(iprintv.ne.0)call tabulate(rss_glm,rss_mod,best_mod,ne_used); call graph(ne_used rss_glm :grid :heading 'Loss in RSS in-sample due to model reduction' :plottype xyplot :nocontact :nolabel :pgborder); call graph(ne_used rss_mod :grid :heading 'Loss in RSS out-of-sample due to model reduction' :plottype xyplot :nocontact :nolabel :pgborder); endif; return; end; call gml_test; b34srun; == ==GLMSOLVE1 Effect of GLM Elasticity alpha on GAS Data Model b34sexec options ginclude('gas.b34'); b34srun; b34sexec matrix; call loaddata; call load(glmsolve :staging); call echooff; lhs='gasout'; rhs='gasout{1 to 20} gasin{1 to 20}'; iprint=2; igraph=1; alower=.1; aupper=.9; na=30; nlam=400; lam_min=.3; thr =.1e-6 ; ne = 20; holdout=0; na=10; call makeglobal(gasout,gasin); call glmsolve(lhs,rhs,alower,aupper,na,nlam, thr,lam_min,thr,ne,rss,rsq,avalue,lam_val, holdout,iprint,igraph); call makelocal(gasout,gasin); b34srun; == ==GLMSOLVE2 Effect of GLM Elasticity alpha on RES72 Model b34sexec options ginclude('b34sdata.mac') member(res72); b34srun; b34sexec matrix; call loaddata; call load(glmsolve :staging); call echooff; lhs='lnq'; rhs='lnl lnk time lnrm2'; iprint=2; igraph=1; alower=.1; aupper=.9; na=30; nlam=400; lam_min=.3; thr =.1e-6 ; ne = 20; holdout=0; na=10; call makeglobal(lnq lnl lnk time lnrm2); call glmsolve(lhs,rhs,alower,aupper,na,nlam, thr,lam_min,thr,ne,rss,rsq,avalue,lam_val,holdout, iprint,igraph); call makelocal(lnq lnl lnk time lnrm2); b34srun; == ==GLS Nonlinear Least Squares applied to GLS Model %b34slet dob34s1=1; %b34slet dob34s2=1; %b34slet dorats =0; %b34slet dosas1 =0; %b34slet dosas2 =0; /$ /$ ********************************************** /$ For AR(1) Models /$ SAS ML like RATS ML /$ B34S GLS like Rats CORC /$ SAS ULS quite different rho(1) /$ ********************************************** /$ b34sexec options ginclude('gas.b34'); b34srun; %b34sif(&dob34s1.eq.1)%then; b34sexec regression maxgls=3; model gasout=gasin; b34srun; b34sexec matrix; call loaddata; call load(gls :staging); call print(gls,glsmod); nn=6; call olsq(gasout gasin{1 to nn} gasout{1 to nn} :print :savex); call glsset; %maxgls=3; %nl2sol=1; %plot=1; call gls; b34srun; %b34sendif; %b34sif(&dosas1.eq.1)%then; /$ B34SEXEC OPTIONS OPEN('testsas.sas') UNIT(29) DISP=UNKNOWN$ B34SRUN$ B34SEXEC OPTIONS CLEAN(29) $ B34SEEND$ B34SEXEC PGMCALL IDATA=29 ICNTRL=29$ SAS $ PGMCARDS$ data new; set b34sdata; gasin_1=lag1(gasin); gasin_2=lag2(gasin); gasin_3=lag3(gasin); gasin_4=lag4(gasin); gasin_5=lag5(gasin); gasin_6=lag6(gasin); gasout_1=lag1(gasout); gasout_2=lag2(gasout); gasout_3=lag3(gasout); gasout_4=lag4(gasout); gasout_5=lag5(gasout); gasout_6=lag6(gasout); run; proc autoreg; model gasout = gasin_1 gasin_2 gasin_3 gasin_4 gasin_5 gasin_6 gasout_1 gasout_2 gasout_3 gasout_4 gasout_5 gasout_6 / method=uls nlag=3; model gasout = gasin_1 gasin_2 gasin_3 gasin_4 gasin_5 gasin_6 gasout_1 gasout_2 gasout_3 gasout_4 gasout_5 gasout_6 / method=ml nlag=3; model gasout = gasin_1 gasin_2 gasin_3 gasin_4 gasin_5 gasin_6 gasout_1 gasout_2 gasout_3 gasout_4 gasout_5 gasout_6 / method=yw nlag=3; model gasout = gasin_1 gasin_2 gasin_3 gasin_4 gasin_5 gasin_6 gasout_1 gasout_2 gasout_3 gasout_4 gasout_5 gasout_6 / method=ityw nlag=3 maxiter=900; run; B34SRETURN$ B34SRUN $ B34SEXEC OPTIONS CLOSE(29)$ B34SRUN$ /$ The next card has to be modified to point to SAS location /$ Be sure and wait until SAS gets done before letting B34S resume B34SEXEC OPTIONS dodos('start /w /r sas testsas') dounix('sas testsas')$ B34SRUN$ B34SEXEC OPTIONS NPAGEOUT NOHEADER WRITEOUT(' ','Output from SAS',' ',' ') WRITELOG(' ','Output from SAS',' ',' ') COPYFOUT('testsas.lst') COPYFLOG('testsas.log') dodos('erase testsas.sas','erase testsas.lst','erase testsas.log') dounix('rm testsas.sas','rm testsas.lst','rm testsas.log')$ B34SRUN$ B34SEXEC OPTIONS HEADER$ B34SRUN$ %b34sendif; /$ /$ user must change parameters market with => /$ b34sexec options ginclude('b34sdata.mac') macro(res72)$ b34srun$ %b34sif(&dob34s2.eq.1)%then; b34sexec regression maxgls=1; model lnq=lnk lnl lnrm2 time; b34srun; b34sexec matrix; call loaddata; call load(gls :staging); call olsq(lnq lnk lnl lnrm2 time :print :savex); call glsset; %maxgls=1; %nl2sol=1; %plot=1; call gls; call olsq(lnq lnk lnl lnrm2 time :print :savex); call glsset; %maxgls=2; %nl2sol=0; %plot=1; call gls; b34srun; %b34sendif; %b34sif(&dorats.eq.1)%then; B34SEXEC OPTIONS OPEN('rats.dat') UNIT(28) DISP=UNKNOWN$ B34SRUN$ B34SEXEC OPTIONS OPEN('rats.in') UNIT(29) DISP=UNKNOWN$ B34SRUN$ B34SEXEC OPTIONS CLEAN(28)$ B34SRUN$ B34SEXEC OPTIONS CLEAN(29)$ B34SRUN$ B34SEXEC PGMCALL$ RATS PASSASTS PCOMMENTS('* ', '* Data passed from B34S(r) system to RATS', '* ') $ PGMCARDS$ * * Various options of RATS AR1 command are tested * linreg lnq # constant lnk lnl lnrm2 time * ar1(method=corc) lnq # constant lnk lnl lnrm2 time * ar1(method=hilu) lnq # constant lnk lnl lnrm2 time * ar1(method=maxl) lnq # constant lnk lnl lnrm2 time * ar1(method=search) lnq # constant lnk lnl lnrm2 time * B34SRETURN$ B34SRUN $ B34SEXEC OPTIONS CLOSE(28)$ B34SRUN$ B34SEXEC OPTIONS CLOSE(29)$ B34SRUN$ B34SEXEC OPTIONS /$ dodos('start /w /r rats386 rats.in rats.out ') dodos('start /w /r rats32s rats.in /run') dounix('rats rats.in rats.out')$ B34SRUN$ B34SEXEC OPTIONS NPAGEOUT WRITEOUT('Output from RATS',' ',' ') COPYFOUT('rats.out') dodos('ERASE rats.in','ERASE rats.out','ERASE rats.dat') dounix('rm rats.in','rm rats.out','rm rats.dat') $ B34SRUN$ %b34sendif; %b34sif(&dosas2.eq.1)%then; /$ B34SEXEC OPTIONS OPEN('testsas.sas') UNIT(29) DISP=UNKNOWN$ B34SRUN$ B34SEXEC OPTIONS CLEAN(29) $ B34SEEND$ B34SEXEC PGMCALL IDATA=29 ICNTRL=29$ SAS $ PGMCARDS$ proc autoreg; model lnq = lnk lnl lnrm2 time / method=uls nlag=1; proc autoreg; model lnq = lnk lnl lnrm2 time / method=ml nlag=1; proc autoreg; model lnq = lnk lnl lnrm2 time / method=yw nlag=1; proc autoreg; model lnq = lnk lnl lnrm2 time / method=ityw maxiter =900 nlag=1; proc autoreg; model lnq = lnk lnl lnrm2 time / method=uls nlag=2; proc autoreg; model lnq = lnk lnl lnrm2 time / method=ml nlag=2; proc autoreg; model lnq = lnk lnl lnrm2 time / method=yw nlag=2; proc autoreg; model lnq = lnk lnl lnrm2 time / method=ityw maxiter =900 nlag=2; run; B34SRETURN$ B34SRUN $ B34SEXEC OPTIONS CLOSE(29)$ B34SRUN$ /$ The next card has to be modified to point to SAS location /$ Be sure and wait until SAS gets done before letting B34S resume B34SEXEC OPTIONS dodos('start /w /r sas testsas') dounix('sas testsas')$ B34SRUN$ B34SEXEC OPTIONS NPAGEOUT NOHEADER WRITEOUT(' ','Output from SAS',' ',' ') WRITELOG(' ','Output from SAS',' ',' ') COPYFOUT('testsas.lst') COPYFLOG('testsas.log') dodos('erase testsas.sas','erase testsas.lst','erase testsas.log') dounix('rm testsas.sas','rm testsas.lst','rm testsas.log')$ B34SRUN$ B34SEXEC OPTIONS HEADER$ B34SRUN$ %b34sendif; == ==GLESJER Glesjer (1969) Heteroscedasticity %b34slet runrats=0; b34sexec options ginclude('greene.mac') member(a5_1); b34srun; b34sexec matrix; call loaddata; call load(g_quandt :staging); call load(b_g_test ); call load(het_test :staging); call load(wald :staging); call load(glesjer :staging); call echooff; mmm =(exp .gt. 0.0); call olsq(exp age ownrent income incomesq :sample mmm :print :savex); /; Glesjer (1969) test See Greene(4th) Page 511 iprint1=0; iprint2=1; /; As per Greene we take only income and incomesq and run the test. /; Note that %x has the sub sample data. xdata=catcol(%x(,3),%x(,4)); call glesjer(%res,xdata,glesjer,prob,typetest,iprint1,iprint2); b34srun; %b34sif(&runrats.eq.1)%then; b34sexec data set; rename exp=exp2; b34srun; b34sexec options open('rats.dat') unit(28) disp=unknown$ b34srun$ b34sexec options open('rats.in') unit(29) disp=unknown$ b34srun$ b34sexec options clean(28)$ b34srun$ b34sexec options clean(29)$ b34srun$ b34sexec pgmcall$ rats passasts pcomments('* ', '* Data passed from B34S(r) system to RATS', '* ', "display @1 %dateandtime() @33 ' Rats Version ' %ratsversion()" '* ') $ PGMCARDS$ * * Rats verification of Glesjer test in B34S * set exp_gt_0 = exp2>0.0 smpl(series=exp_gt_0) linreg exp2 / res1 # constant age ownrent income incomesq set abs_res1 = abs(res1) set res1_res1 = res1*res1 set log_abs_res1 = log(abs_res1) * print * * exp2 age ownrent income incomesq res1 abs_res1 $ * res1_res1 log_abs_res1 linreg res1_res1 # constant income incomesq test # 2 3 # 0.0 0.0 linreg(robust) res1_res1 # constant income incomesq test # 2 3 # 0.0 0.0 linreg abs_res1 # constant income incomesq test # 2 3 # 0.0 0.0 linreg(robust) abs_res1 # constant income incomesq test # 2 3 # 0.0 0.0 linreg log_abs_res1 # constant income incomesq test # 2 3 # 0.0 0.0 linreg(robust) log_abs_res1 # constant income incomesq test # 2 3 # 0.0 0.0 b34sreturn$ b34srun $ b34sexec options close(28)$ b34srun$ b34sexec options close(29)$ b34srun$ b34sexec options /$ dodos(' rats386 rats.in rats.out ') dodos('start /w /r rats32s rats.in /run') dounix('rats rats.in rats.out')$ B34SRUN$ b34sexec options npageout WRITEOUT('Output from RATS',' ',' ') COPYFOUT('rats.out') dodos('ERASE rats.in','ERASE rats.out','ERASE rats.dat') dounix('rm rats.in','rm rats.out','rm rats.dat') $ B34SRUN$ == ==GLESJER2 Glesjer (1969) Heteroscedasticity /; /; Example From Greene (Edition 4) page 510 /; Six tests. Rats varifed results /; %b34slet runrats=0; b34sexec options ginclude('greene.mac') member(a5_1); b34srun; b34sexec matrix; call loaddata; call load(g_quandt :staging); call load(b_g_test ); call load(het_test :staging); call load(wald :staging); call load(glesjer :staging); call echooff; mmm =(exp .gt. 0.0); /; Parameters for het_test call olsq(exp age ownrent income incomesq :sample mmm :print :qr :savex); /; Parameters for het_test res1=%res; xhold=%x; bg_order=4; n1=36; n2=37; gqprint1=0; gqprint2=1; isort1=3; isort2=3; call het_test(%res,%y,%x,%names,%lag,bg_order,1,1,0,1,1,0,isort1,isort2, n1,n2,gqprint1,gqprint2); call olsq(exp income incomesq :sample mmm :print :savex); reshold=%res; hxpxinv=%xpxinv; xhold=%x; isort1=1; isort2=1; call het_test(%res,%y,%x,%names,%lag,bg_order,1,1,0,1,1,0,isort1,isort2, n1,n2,gqprint1,gqprint2); /; Glesjer (1969) test See Greene(4th) Page 511 iprint1=0; iprint2=1; res=res1; xdata=catcol(xhold(,1),xhold(,2)); call glesjer(res,xdata,glesjer,prob,typetest,iprint1,iprint2); b34srun; %b34sif(&runrats.eq.1)%then; b34sexec data set; rename exp=exp2; b34srun; b34sexec options open('rats.dat') unit(28) disp=unknown$ b34srun$ b34sexec options open('rats.in') unit(29) disp=unknown$ b34srun$ b34sexec options clean(28)$ b34srun$ b34sexec options clean(29)$ b34srun$ b34sexec pgmcall$ rats passasts pcomments('* ', '* Data passed from B34S(r) system to RATS', '* ', "display @1 %dateandtime() @33 ' Rats Version ' %ratsversion()" '* ') $ PGMCARDS$ * * Rats verification of Glesjer test in B34S * set exp_gt_0 = exp2>0.0 smpl(series=exp_gt_0) linreg exp2 / res1 # constant age ownrent income incomesq set abs_res1 = abs(res1) set res1_res1 = res1*res1 set log_abs_res1 = log(abs_res1) * print * * exp2 age ownrent income incomesq res1 abs_res1 $ * res1_res1 log_abs_res1 linreg res1_res1 # constant income incomesq test # 2 3 # 0.0 0.0 linreg(robust) res1_res1 # constant income incomesq test # 2 3 # 0.0 0.0 linreg abs_res1 # constant income incomesq test # 2 3 # 0.0 0.0 linreg(robust) abs_res1 # constant income incomesq test # 2 3 # 0.0 0.0 linreg log_abs_res1 # constant income incomesq test # 2 3 # 0.0 0.0 linreg(robust) log_abs_res1 # constant income incomesq test # 2 3 # 0.0 0.0 b34sreturn$ b34srun $ b34sexec options close(28)$ b34srun$ b34sexec options close(29)$ b34srun$ b34sexec options /$ dodos(' rats386 rats.in rats.out ') dodos('start /w /r rats32s rats.in /run') dounix('rats rats.in rats.out')$ B34SRUN$ b34sexec options npageout WRITEOUT('Output from RATS',' ',' ') COPYFOUT('rats.out') dodos('ERASE rats.in','ERASE rats.out','ERASE rats.dat') dounix('rm rats.in','rm rats.out','rm rats.dat') $ B34SRUN$ == ==GLSDATA GLS b34sexec options ginclude('b34sdata.mac') macro(res72)$ b34srun$ b34sexec regression maxgls=1 toll=.1e-14; model lnq=lnk lnl lnrm1 time; b34srun; b34sexec matrix; call loaddata; call load(glsdata :staging); call echooff; /$ subroutine glsdata(rho,x,newx); /$ /$ Transforms x or y given rho the GLS parameter /$ X is the data where each col is a series /$ newx is the transformed data /$ /$ GLS using "ols method" call olsq(lnq lnk lnl lnrm1 time :print :savex); acf1=acf(%res,4); adj=dfloat(norows(%res))/ dfloat(integers(norows(%res)-1,norows(%res)-norows(acf1),-1)); acf2=afam(acf1)*afam(adj); call print('acf large sample ',acf(%res,4), 'acf_small ',acf2); %basex=%x; %basey=%y; %baseres=%res; %holdnam=%names; %holdlags=%lag; call olsq(%baseres %baseres{1} :noint :print); call glsdata(%coef,%basey,newy); call glsdata(%coef,%basex,newx); call olsq(newy newx :noint :print); call tabulate(%holdnam %holdlags %coef %se %t); b34srun; == ==GLS_ML GLS AR(1) Estimation using ML %b34slet dosas =0; %b34slet dorats=0; %b34slet limdep=0; %b34slet do_GLS_ML =1; %b34slet do_GLS_co =1; %b34slet do_GLS_2P =1; b34sexec options ginclude('greene.mac') macro(nf2_2)$ b34srun$ b34sexec data set$ build log_gas log_pg log_y log_pnc log_puc log_ppt log_pd log_pn log_ps gasnew log_gnew ; rename g=gas; gen year=year-1959$ /$gen y =y/8513.52 $ gen log_gas =log(gas)$ gen log_pg =log(pg)$ gen log_y =log(y)$ gen log_pnc =log(pnc)$ gen log_puc =log(puc)$ gen log_ppt =log(ppt)$ gen log_pd =log(pd)$ gen log_pn =log(pn)$ gen log_ps =log(ps)$ gen gasnew = gas/pop; label gasnew = 'gasnew = gas/pop'; gen log_gnew = log(gasnew); label log_gnew = ' log_gnew=log(gasnew)'; label log_gas ='log(gas) '$ label log_pg ='log(pg) '$ label log_y ='log(y) '$ label log_pnc ='log(pnc) '$ label log_puc ='log(puc) '$ label log_ppt ='log(ppt) '$ label log_pd ='log(pd) '$ label log_pn ='log(pn) '$ label log_ps ='log(ps) '$ b34seend$ %b34sif(&do_gls_2P.ne.0)%then; b34sexec regression maxgls=1; model log_gnew=log_pg log_y log_pnc log_puc; b34srun; %b34sendif; b34sexec matrix; call loaddata; call load(gls :staging); call load(gls_ml :staging); /$ ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /$ This does Cochrane-Orcutt Method for AR(1) and AR(2) /$ %b34sif(&do_gls_co.ne.0)%then; call olsq(log_gnew log_pg log_y log_pnc log_puc :print :savex); call glsset; %maxgls=1; %plot=1; %nl2sol=1; call gls; call olsq(log_gnew log_pg log_y log_pnc log_puc :print :savex); call glsset; %maxgls=2; %plot=1; %nl2sol=1; call gls; %b34sendif; /$ ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /$ This does Maximum Likehood estimation of AR(1) Model /$ %b34sif(&do_gls_ml.ne.0)%then; call olsq(log_gnew log_pg log_y log_pnc log_puc :print :savex); call glsmlset; call gls_ml; %b34sendif; b34srun; %b34sif(&dosas.eq.1)%then; /$ B34SEXEC OPTIONS OPEN('testsas.sas') UNIT(29) DISP=UNKNOWN$ B34SRUN$ B34SEXEC OPTIONS CLEAN(29) $ B34SEEND$ B34SEXEC PGMCALL IDATA=29 ICNTRL=29$ SAS $ PGMCARDS$ proc autoreg; model log_gnew=log_pg log_y log_pnc log_puc / method=uls nlag=1; model log_gnew=log_pg log_y log_pnc log_puc / method=ml nlag=1; model log_gnew=log_pg log_y log_pnc log_puc / method=yw nlag=1; model log_gnew=log_pg log_y log_pnc log_puc / method=ityw nlag=1 maxiter=900; run; B34SRETURN$ B34SRUN $ B34SEXEC OPTIONS CLOSE(29)$ B34SRUN$ /$ The next card has to be modified to point to SAS location /$ Be sure and wait until SAS gets done before letting B34S resume B34SEXEC OPTIONS dodos('start /w /r sas testsas') dounix('sas testsas')$ B34SRUN$ B34SEXEC OPTIONS NPAGEOUT NOHEADER WRITEOUT(' ','Output from SAS',' ',' ') WRITELOG(' ','Output from SAS',' ',' ') COPYFOUT('testsas.lst') COPYFLOG('testsas.log') dodos('erase testsas.sas','erase testsas.lst','erase testsas.log') dounix('rm testsas.sas','rm testsas.lst','rm testsas.log')$ B34SRUN$ B34SEXEC OPTIONS HEADER$ B34SRUN$ %b34sendif; %b34sif(&dorats.eq.1)%then; B34SEXEC OPTIONS OPEN('rats.dat') UNIT(28) DISP=UNKNOWN$ B34SRUN$ B34SEXEC OPTIONS OPEN('rats.in') UNIT(29) DISP=UNKNOWN$ B34SRUN$ B34SEXEC OPTIONS CLEAN(28)$ B34SRUN$ B34SEXEC OPTIONS CLEAN(29)$ B34SRUN$ B34SEXEC PGMCALL$ RATS PASSASTS PCOMMENTS('* ', '* Data passed from B34S(r) system to RATS', '* ') $ PGMCARDS$ * * Various options of RATS AR1 command are tested * linreg log_gnew # constant log_pg log_y log_pnc log_puc * ar1(method=corc) log_gnew # constant log_pg log_y log_pnc log_puc * ar1(method=hilu) log_gnew # constant log_pg log_y log_pnc log_puc * ar1(method=maxl) log_gnew # constant log_pg log_y log_pnc log_puc * ar1(method=search) log_gnew # constant log_pg log_y log_pnc log_puc * B34SRETURN$ B34SRUN $ B34SEXEC OPTIONS CLOSE(28)$ B34SRUN$ B34SEXEC OPTIONS CLOSE(29)$ B34SRUN$ B34SEXEC OPTIONS /$ dodos('start /w /r rats386 rats.in rats.out ') dodos('start /w /r rats32s rats.in /run') dounix('rats rats.in rats.out')$ B34SRUN$ B34SEXEC OPTIONS NPAGEOUT WRITEOUT('Output from RATS',' ',' ') COPYFOUT('rats.out') dodos('ERASE rats.in','ERASE rats.out','ERASE rats.dat') dounix('rm rats.in','rm rats.out','rm rats.dat') $ B34SRUN$ %b34sendif; %b34sif(&limdep.eq.1)%then; b34sexec options open('limdep.in') unit(29) disp=unknown$ b34SRUN$ b34sexec options clean(29)$ b34srun$ b34sexec pgmcall$ limdep ibat necho title('Results from limdep') pcomments('open ; output = limdep.out $') $ /$ /$ LIMDEP commands such as /$ /$ crmo; lhs=gasout; rhs=one,gasin$ /$ armax; lhs = gasout; rhs= gasin ; model = 3,1,1 $ /$ /$ listed after prmcards$ sentence /$ pgmcards$ crmo; lhs=log_gnew;rhs=one,log_pg,log_y,log_pnc,log_puc $ regress;lhs=log_gnew;rhs=one,log_pg,log_y,log_pnc,log_puc;ar1$ regress;lhs=log_gnew;rhs=one,log_pg,log_y,log_pnc,log_puc;ar1;alg=corc$ regress;lhs=log_gnew;rhs=one,log_pg,log_y,log_pnc,log_puc;ar1;alg=MLE$ b34sreturn$ b34srun $ b34sexec options close(29)$ b34srun$ b34sexec options /$ dodos('start /w /r limdep95 < limdep.in /Wl,p195') dodos('limdep limdep.in') dounix('limdep < limdep.in')$ b34srun$ /$system('d:\limdep\blimdep < limdep.in')$ b34srun$ /$system('blimdep < limdep.in')$ b34srun$ b34sexec options npageout writeout('Output from LIMDEP',' ',' ') copyfout('limdep.out') writelog('Output from LIMDEP log') copyflog('trace.lim') dodos('erase limdep.in ', 'erase limdep.out', 'erase trace.lim ', 'cls') dounix('rm limdep.in ', 'rm limdep.out', 'rm trace.lim '); b34srun$ == ==GLS_2P Two Pass GLS /; GLS Experiments Testing OLS GLS Model Estimation b34sexec options ginclude('b34sdata.mac') macro(res72)$ b34srun$ /; Model B option using REGRESSION Command /; b34sexec regression maxgls=2 toll=.1e-14 ntac=7; /; model lnq=lnk lnl lnrm1 time; b34srun; b34sexec matrix; call loaddata; call load(glsdata :staging); call load(gls_2p :staging); call load(data_acf); call echooff; call olsq(lnq lnk lnl lnrm1 time :print :savex); call gls2pset; %maxgls=2; %doplots=1; %doacf =12; call gls_2p; b34srun; == ==GPH Geweke & Porter-Hudak (1983) Fractional Diff. Test /$ /$ Test for GPH( ) /$ %b34slet runrats =0; b34sexec options ginclude('greene.mac') member(nex18_2); b34srun; /; /; ACF agrees with Greene (2003, 648) but periodogram does not unless /; scaled. Rats code does not agree with Greene but does with B34S /; Greene in e-mail to H. H. Stokes 5 September 2006 agrees with the /; Rats/B34S answer of |.24430| /; b34sexec matrix; call loaddata; call load(gph :staging); call echooff; series1 = dif(log(realgnp)); /; series2 = dif(log(realgnpg)); /; acf1=acf(series1 ,12); /; acf2=acf(series2 ,12); /; call print('Attempting to Replicate Greene (2000, 787)':); /; call tabulate(acf1,acf2); call gph(series1,.51,d,se,se2,1); /; call gph(series1,.51,d,se,se2,2); b34srun; %b34sif(&runrats.ne.0)%then; B34SEXEC OPTIONS OPEN('rats.dat') UNIT(28) DISP=UNKNOWN$ B34SRUN$ B34SEXEC OPTIONS OPEN('rats.in') UNIT(29) DISP=UNKNOWN$ B34SRUN$ B34SEXEC OPTIONS CLEAN(28)$ B34SRUN$ B34SEXEC OPTIONS CLEAN(29)$ B34SRUN$ B34SEXEC PGMCALL$ RATS PASSASTS pcomments('* ', '* Data passed from B34S(r) system to RATS', '* ', "display @1 %dateandtime() @33 ' Rats Version ' %ratsversion()" '* ') $ PGMCARDS$ * SOURCE(NOECHO) d:\R\GPH.SRC * SOURCE(ECHO) d:\R\GPH.SRC * set y = log(realgnp) SET SERIES = y - y{1} * Frequencies used = T**POWER. Default = .5 * .51 replicated greene's # of freq * For detail see Geweke-Porter-Hudak @gph(power=.51) series B34SRETURN$ B34SRUN $ B34SEXEC OPTIONS CLOSE(28)$ B34SRUN$ B34SEXEC OPTIONS CLOSE(29)$ B34SRUN$ B34SEXEC OPTIONS /$ dodos('start /w /r rats386 rats.in rats.out ') dodos('start /w /r rats32s rats.in /run') dounix('rats rats.in rats.out')$ B34SRUN$ B34SEXEC OPTIONS NPAGEOUT WRITEOUT('Output from RATS',' ',' ') COPYFOUT('rats.out') dodos('ERASE rats.in','ERASE rats.out','ERASE rats.dat') dounix('rm rats.in','rm rats.out','rm rats.dat') $ B34SRUN$ %b34sendif; == ==G_QUANDT Goldfeld-Quandt Heteroskedasticity Test /; /; Example From Greene (Edition 4) page 510 /; Greene (Edition 5) page 224 /; /; Greene answer is 15.001 /; b34sexec options ginclude('greene.mac') member(a5_1); b34srun; b34sexec matrix; call loaddata; call load(g_quandt :staging); call echooff; call print('Example From Greene (Edition 4) page 510':); call print(' Greene (Edition 5) page 224':); call print('Greene answer is 15.001':); mmm =(exp .gt. 0.0); call olsq(exp age ownrent income incomesq :sample mmm :print :savex); /; call print(%y,%x); /; Help on arguments to test /; Goldfield-Quandt Heteroscedasticity Test /; Goldfield-Quandt "Some Tests for Homoscedasticity," /; Journal of the American Statistical Association /; 60, 1965, pp539-547 /; /; y => left hand variable /; x => Right hand variable matrix /; sortcol => If ne 0 sets sort col /; names => Names Col /; lag => Lag Col /; n1 => Last obs in left hand side /; n2 => First Obs in right hand side /; ss1 => e'e for left hand side /; ss2 => e'e for right hand side /; g_q_test => Goldfield-Quandt test /; prob => Probability of test /; iprint1 => Print regressions /; iprint2 => Print test do i=3,3; call g_quandt(%y,%x,i,%names,%lag,36,37,ss1,ss2,g_q,prob,1,1); enddo; b34srun; == ==G_QUANDT2 Sorts over ranges and with alt # dropped. b34sexec options ginclude('b34sdata.mac') member(grunfeld); b34srun; b34sexec matrix; call loaddata; call load(g_quandt :staging); call echooff; /; call print(%y,%x); /; Help on arguments to test /; Goldfield-Quandt Heteroscedasticity Test /; Goldfield-Quandt "Some Tests for Homoscedasticity," /; Journal of the American Statistical Association /; 60, 1965, pp539-547 /; /; y => left hand variable /; x => Right hand variable matrix /; sortcol => If ne 0 sets sort col /; names => Names Col /; lag => Lag Col /; n1 => Last obs in left hand side /; n2 => First Obs in right hand side /; ss1 => e'e for left hand side /; ss2 => e'e for right hand side /; g_q_test => Goldfield-Quandt test /; prob => Probability of test /; iprint1 => Print regressions /; iprint2 => Print test call print(' ':); call print('GE Equation':); call print('-----------':); call print(' ':); call olsq(i_ge f_ge c_ge :diag :print :savex); call print(' ':); call describe(%res :print); /; call print(%x); nn1=10; nn2=11; do i=1,2; do j=0,4; n1=nn1-j; n2=nn2+j; call g_quandt(%y,%x,i,%names,%lag,n1,n2,ss1,ss2,g_q,prob,1,1); enddo; enddo; b34srun; == ==HETTEST1 Theil Heteroskedasticity Test on Sorted Data b34sexec options ginclude('greene08.mac') member(table25_1); b34srun; b34sexec regression residuala; model majordrg=age income exp_inc avgexp ownrent selfempl depndt inc_per cur_add major active; b34srun; b34sexec sort; by age; b34srun; b34sexec regression residuala; model majordrg=age income exp_inc avgexp ownrent selfempl depndt inc_per cur_add major active; b34srun; /; b34sexec matrix; call loaddata; call load(hettest1 :staging); call echooff; call olsq(majordrg age income exp_inc avgexp ownrent selfempl depndt inc_per cur_add major active :print :savex); call print(' ':); isort1=0; isort2=0; pertest=.33333333; print1=1; call hettest1(%res,%x,%names,%lag,f,rf,sig,isort1,isort2, pertest,print1); isort1=1; isort2=9; pertest=.33333333; print1=1; call hettest1(%res,%x,%names,%lag,f,rf,sig,isort1,isort2, pertest,print1); b34srun; == ==HET_TEST1_2 Theil & Goldfield Quandt tests b34sexec options ginclude('greene08.mac') member(table25_1); b34srun; b34sexec regression residuala; model majordrg=age income exp_inc avgexp ownrent selfempl depndt inc_per cur_add major active; b34srun; b34sexec sort; by age; b34srun; b34sexec regression residuala; model majordrg=age income exp_inc avgexp ownrent selfempl depndt inc_per cur_add major active; b34srun; /; b34sexec matrix; call loaddata; call load(hettest1 :staging); call load(g_quandt :staging); call echooff; call olsq(majordrg age income exp_inc avgexp ownrent selfempl depndt inc_per cur_add major active :print :savex); call print(' ':); isort1=0; isort2=0; pertest=.33333333; print1=1; call hettest1(%res,%x,%names,%lag,f,rf,sig,isort1,isort2, pertest,print1); isort1=1; isort2=4; pertest=.33333333; print1=1; /; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ call hettest1(%res,%x,%names,%lag,f,rf,sig,isort1,isort2, pertest,print1); /; Help on arguments to test /; Goldfield-Quandt Heteroscedasticity Test /; Goldfield-Quandt "Some Tests for Homoscedasticity," /; Journal of the American Statistical Association /; 60, 1965, pp539-547 /; /; y => left hand variable /; x => Right hand variable matrix /; sortcol => If ne 0 sets sort col /; names => Names Col /; lag => Lag Col /; n1 => Last obs in left hand side /; n2 => First Obs in right hand side /; ss1 => e'e for left hand side /; ss2 => e'e for right hand side /; g_q_test => Goldfield-Quandt test /; prob => Probability of test /; iprint1 => Print regressions /; iprint2 => Print test do i=1,4; n1=norows(%y)/3; n2=norows(%y)-n1; call g_quandt(%y,%x,i,%names,%lag,n1,n2,ss1,ss2,g_q,prob,0,1); enddo; b34srun; == ==HET_TEST Heteroskedasticity Tests /; /; Example From Greene (Edition 4) page 510 /; Example From Greene (Edition 5) page 224-225 /; b34sexec options ginclude('greene.mac') member(a5_1); b34srun; b34sexec matrix; call loaddata; call load(g_quandt :staging); call load(b_g_test ); call load(het_test :staging); call echooff; mmm =(exp .gt. 0.0); call olsq(exp age ownrent income incomesq :sample mmm :print :savex); resfull=%res; /; call print(%y,%x); /; Parameters for het_test n1=36; n2=37; gqprint1=0; gqprint2=1; isort1=3; isort2=3; /; subroutine het_test(res,y,x,names,lag,bp1,bp2,bp3, /; white1,white2,white3, /; isort1,isort2,n1,n2,gqprint1,gqprint2); /; /; res => Residual vector usually %res /; y => left hand side usually %y /; x => x matrix %x /; names => Used for Goldfeld-Quandt test /; lags => Used for Goldfeld-Quandt test /; bp1 => 1 do t variant of Breusch Pagan Test /; bp2 => 1 do x variant of Breusch-Pagan Test /; bp3 => 1 do x**2 variant of Breusch-Pagan Test /; white1 => 1 do t variant of White Test /; white2 => 1 do x variant of White Test /; white3 => 1 do x**2 variant of White Test /; isort1 => Beginning col to sort for Goldfeld-Quandt test /; isort2 => End col to sort for Goldfeld-Quandt test /; if isort1 or isort2 = 0 => no sorting is done. /; n1 => last obs for top of Goldfeld-Quandt test /; n2 => first obs for top of Goldfeld-Quandt test /; gqprint1 => print regressions for Goldfled-Quandt test /; gqprint2 => print results of Goldfeld-Quandt test /; /; Breusch-Pagan (1979) Test Three variants - Hetroscedasticity /; White Test Three variants - Hetroscedasticity /; Goldfeld-Quandt(1965) Requires sort index /; /; Routine built December 2005. Arguments subject to change /; This gets X call print('Tests run on all variables on right':); call print('This gets the Greene Goldfield Quandt value':); call print('+++++++++++++++++++++++++++++++++++++++++++':); call het_test(resfull,%y,%x,%names,%lag,1,1,0,1,1,0,isort1,isort2, n1,n2,gqprint1,gqprint2); call print('These tests restrict X as per Greene ':); call print('Will obtain Greene Breusch-Pagan Value':); call print('Note effect on Goldfield-Quandt and Koenker-Bassett':); call print('+++++++++++++++++++++++++++++++++++++':); call olsq(exp income incomesq :sample mmm :print :savex); isort1=1; isort2=1; call het_test(resfull,%y,%x,%names,%lag,1,1,0,1,1,0,isort1,isort2, n1,n2,gqprint1,gqprint2); b34srun; == ==IACF Inverse ACF %b34sif(&dorats.ne.0)%then; b34sexec options open('rats.dat') unit(28) disp=unknown$ b34srun$ b34sexec options open('rats.in') unit(29) disp=unknown$ b34srun$ b34sexec options clean(28)$ b34srun$ b34sexec options clean(29)$ b34srun$ b34sexec pgmcall$ rats passasts pcomments('* ', '* Data passed from B34S(r) system to RATS', '* ', "display @1 %dateandtime() @33 ' Rats Version ' %ratsversion()" '* ') $ PGMCARDS$ * correlate(inverse=iacf,partial=pacf) gasin b34sreturn$ b34srun $ b34sexec options close(28)$ b34srun$ b34sexec options close(29)$ b34srun$ b34sexec options /$ dodos(' rats386 rats.in rats.out ') dodos('start /w /r rats32s rats.in /run') dounix('rats rats.in rats.out')$ B34SRUN$ b34sexec options npageout WRITEOUT('Output from RATS',' ',' ') COPYFOUT('rats.out') dodos('ERASE rats.in','ERASE rats.out','ERASE rats.dat') dounix('rm rats.in','rm rats.out','rm rats.dat') $ B34SRUN$ %b34sendif; == ==IACF1 Test data from Wei /; Note widely different answers /; Four programs being tested /; iacf => Uses Wei setup and estimated the AR model with OLS /; Wei(2006) "Time Series Analysis" /; /; iacfn => Uses Bovas Abraham & Johannes Ledolter model using PACF /; Bovas Abraham & Johannes Ledolter /; "A Note on Inverse Autocorrelations" /; Biometrika (1984) Vol 71 # 3 pp 609-614 eq 1-5 /; /; iacf(k) = (-ar(k) + sum (j=1,p-k) ar(j)*ar(j+k)) / /; (1.0 +sumsq(ar)) /; /; where ar = partial acf /; /; %b34slet dosca=1; %b34slet dosas=0; %b34slet dorats=1; b34sexec options ginclude('wei.mac') member(wei_w1); b34srun; b34sexec matrix; call loaddata; call load(iacf :staging); call load(iacfn :staging); call echooff; call print('+++++++++++++++++++++++++++++':); call print('Test Case From Wei page 127 ':); call print('Note that SCA and B34S get same answers but Wei differs':); n=10; acf1= acf(defects,n,se,pacf); call iacf(defects,iacf1,n); call iacfn(defects,iacf2,n); call tabulate(acf1,pacf,iacf1,iacf2 :title 'ACF, PACF IACF iacfn'); b34srun; %b34sif(&dosca.ne.0)%then; b34sexec options open('sca.dat') disp=unknown unit(28)$ b34srun$ b34sexec options open('sca.cmd') disp=unknown unit(29)$ b34srun$ b34sexec options clean(28)$ b34srun$ b34sexec options clean(29)$ b34srun$ b34sexec pgmcall$ sca scafname=mydata$ pgmcards$ /$#==myrun --- these commands are required to load the b34s data. --- assign file 18. attrib access(read). external 'sca.dat'. call procedure is mydata. file is 18. --- iacf defects stop. return /$#== b34sreturn$ b34srun$ b34sexec options close(28)$ b34srun$ b34sexec options close(29)$ b34srun$ b34sexec options dounix('sca sca.cmd > sca.out') dodos('scaw32 10000 /f:sca.cmd /p:myrun /o:sca.out') $ b34srun$ b34sexec options npageout writeout('output from sca',' ',' ') copyfout('sca.out') dodos('erase sca.cmd','erase sca.out','erase sca.dat') dounix('rm sca.cmd','rm sca.out','rm sca.dat') $ b34srun$ %b34sendif; %b34sif(&dorats.ne.0)%then; b34sexec options open('rats.dat') unit(28) disp=unknown$ b34srun$ b34sexec options open('rats.in') unit(29) disp=unknown$ b34srun$ b34sexec options clean(28)$ b34srun$ b34sexec options clean(29)$ b34srun$ b34sexec pgmcall$ rats passasts pcomments('* ', '* Data passed from B34S(r) system to RATS', '* ', "display @1 %dateandtime() @33 ' Rats Version ' %ratsversion()" '* ') $ PGMCARDS$ * correlate(inverse=iacf,partial=pacf) defects b34sreturn$ b34srun $ b34sexec options close(28)$ b34srun$ b34sexec options close(29)$ b34srun$ b34sexec options /$ dodos(' rats386 rats.in rats.out ') dodos('start /w /r rats32s rats.in /run') dounix('rats rats.in rats.out')$ B34SRUN$ b34sexec options npageout WRITEOUT('Output from RATS',' ',' ') COPYFOUT('rats.out') dodos('ERASE rats.in','ERASE rats.out','ERASE rats.dat') dounix('rm rats.in','rm rats.out','rm rats.dat') $ B34SRUN$ %b34sendif; %b34sif(&dosas.eq.1)%then; b34sexec options open('testsas.sas') unit(29) disp=unknown$ b34srun$ b34sexec options clean(29) $ b34seend$ b34sexec pgmcall idata=29 icntrl=29$ sas $ * sas commands next ; pgmcards$ proc arima; identify var=defects; run; b34sreturn$ b34srun $ b34sexec options close(29)$ b34srun$ b34sexec options dodos('start /w /r sas testsas' ) dounix('sas testsas' ) $ b34srun$ b34sexec options npageout noheader writeout(' ','output from sas',' ',' ') writelog(' ','output from sas',' ',' ') copyfout('testsas.lst') copyflog('testsas.log') dodos('erase testsas.sas','erase testsas.lst','erase testsas.log') dounix('rm testsas.sas','rm testsas.lst','rm testsas.log') $ b34srun$ b34sexec options header$ b34srun$ %b34sendif; b34sexec options ginclude('wei.mac') member(wei_w7); b34srun; b34sexec data set; build ln_lynx; gen ln_lynx = dlog(lynx); b34srun; b34sexec matrix; call loaddata; call load(iacf :staging); call load(iacfn :staging); call echooff; call print('+++++++++++++++++++++++++++++':); call print('Test Case From Wei page 132 ':); call print('Note that SCA and B34S get same answers but Wei differs':); n=10; acf1= acf(ln_lynx,n,se,pacf); call iacf(ln_lynx,iacf1,n); call iacfn(ln_lynx,iacf2,n); call tabulate(acf1,pacf,iacf1,iacf2 :title 'ACF, PACF IACF iacfn'); b34srun; %b34sif(&dosca.ne.0)%then; b34sexec options open('sca.dat') disp=unknown unit(28)$ b34srun$ b34sexec options open('sca.cmd') disp=unknown unit(29)$ b34srun$ b34sexec options clean(28)$ b34srun$ b34sexec options clean(29)$ b34srun$ b34sexec pgmcall$ sca scafname=mydata$ pgmcards$ /$#==myrun --- these commands are required to load the b34s data. --- assign file 18. attrib access(read). external 'sca.dat'. call procedure is mydata. file is 18. --- acf ln_lynx pacf ln_lynx iacf ln_lynx stop. return /$#== b34sreturn$ b34srun$ b34sexec options close(28)$ b34srun$ b34sexec options close(29)$ b34srun$ b34sexec options dounix('sca sca.cmd > sca.out') dodos('scaw32 10000 /f:sca.cmd /p:myrun /o:sca.out') $ b34srun$ b34sexec options npageout writeout('output from sca',' ',' ') copyfout('sca.out') dodos('erase sca.cmd','erase sca.out','erase sca.dat') dounix('rm sca.cmd','rm sca.out','rm sca.dat') $ b34srun$ %b34sendif; %b34sif(&dorats.ne.0)%then; b34sexec options open('rats.dat') unit(28) disp=unknown$ b34srun$ b34sexec options open('rats.in') unit(29) disp=unknown$ b34srun$ b34sexec options clean(28)$ b34srun$ b34sexec options clean(29)$ b34srun$ b34sexec pgmcall$ rats passasts pcomments('* ', '* Data passed from B34S(r) system to RATS', '* ', "display @1 %dateandtime() @33 ' Rats Version ' %ratsversion()" '* ') $ PGMCARDS$ * correlate(inverse=iacf,partial=pacf) ln_lynx b34sreturn$ b34srun $ b34sexec options close(28)$ b34srun$ b34sexec options close(29)$ b34srun$ b34sexec options /$ dodos(' rats386 rats.in rats.out ') dodos('start /w /r rats32s rats.in /run') dounix('rats rats.in rats.out')$ B34SRUN$ b34sexec options npageout WRITEOUT('Output from RATS',' ',' ') COPYFOUT('rats.out') dodos('ERASE rats.in','ERASE rats.out','ERASE rats.dat') dounix('rm rats.in','rm rats.out','rm rats.dat') $ B34SRUN$ %b34sendif; == ==LASSO Lasso Shrinkage /; /; Ridge regression also shown /; /; For further detail see Stokes (200x, Chapter 10) /; /; %b34slet dorats =0; %b34slet doridge=1; b34sexec options ginclude('gas.b34'); b34srun; b34sexec matrix; call loaddata; call load(ridge :staging); call load(lasso :staging); call echooff; k=6; call olsq(gasout gasout{1 to k} gasin{1 to k} :print :savex); %b34sif(&doridge.ne.0)%then; lamda=2.; call ridge(%y,%x,lamda,%coef,%names,%lag,ridge_c,0); call ridge(%y,%x,lamda,%coef,%names,%lag,ridge_c,1); %b34sendif; lamda=10.*sum(%coef)/2.0; call echooff; call lasso(%y,%x,%coef,%lcoef1,%l_t1,lamda,lresid1,1); call lasso(%y,%x,%coef,%lcoef2,%l_t2,lamda,lresid2,3); call tabulate(%names,%lag,%coef,%se,%t,%lcoef1,%l_t1,%lcoef2,%l_t2); b34srun; %b34sif(&dorats.ne.0)then; b34sexec options open('rats.dat') unit(28) disp=unknown$ b34srun$ b34sexec options open('rats.in') unit(29) disp=unknown$ b34srun$ b34sexec options clean(28)$ b34srun$ b34sexec options clean(29)$ b34srun$ /; Uses logic from RATS User's Guide Version 6.1 Page 192 b34sexec pgmcall$ rats passasts pcomments('* ', '* Data passed from B34S(r) system to RATS', '* ', "display @1 %dateandtime() @33 ' Rats Version ' %ratsversion()" '* ') $ PGMCARDS$ * * Non centered results * cmoment # constant gasout{1 to 6} gasin{1 to 6} gasout linreg(cmoment) gasout # constant gasout{1 to 6} gasin{1 to 6} do row=1,13 compute %cmom(row,row)=%cmom(row,row)+2 end do row linreg(cmoment) gasout # constant gasout{1 to 6} gasin{1 to 6} b34sreturn$ b34srun $ b34sexec options close(28)$ b34srun$ b34sexec options close(29)$ b34srun$ b34sexec options /$ dodos(' rats386 rats.in rats.out ') dodos('start /w /r rats32s rats.in /run') dounix('rats rats.in rats.out')$ B34SRUN$ b34sexec options npageout WRITEOUT('Output from RATS',' ',' ') COPYFOUT('rats.out') dodos('ERASE rats.in','ERASE rats.out','ERASE rats.dat') dounix('rm rats.in','rm rats.out','rm rats.dat') $ B34SRUN$ %b34sendif; == ==LASSO2 Generalized Lasso Shrinkage /; /; Ridge regression also shown /; /; For further detail see Stokes (200x, Chapter 10) /; /; %b34slet dorats =0; %b34slet doridge=1; b34sexec options ginclude('gas.b34'); b34srun; b34sexec matrix; call loaddata; call load(ridge :staging); call load(lasso :staging); call load(lasso2 :staging); call echooff; k=6; call olsq(gasout gasout{1 to k} gasin{1 to k} :print :savex); %b34sif(&doridge.ne.0)%then; lamda=2.; call ridge(%y,%x,lamda,%coef,%names,%lag,ridge_c,0); call ridge(%y,%x,lamda,%coef,%names,%lag,ridge_c,1); %b34sendif; lamda=10.*sum(%coef)/2.0; call echooff; call lasso(%y,%x,%coef,%lcoef1,%l_t1,lamda,lresid1,1); call lasso(%y,%x,%coef,%lcoef2,%l_t2,lamda,lresid2,3); call tabulate(%names,%lag,%coef,%se,%t,%lcoef1,%l_t1,%lcoef2,%l_t2); call lasso2(%y,%x,%coef,%lcoef1,%l_t1,lamda,lresid1,1,2.); call lasso2(%y,%x,%coef,%lcoef2,%l_t2,lamda,lresid2,3,2.); call tabulate(%names,%lag,%coef,%se,%t,%lcoef1,%l_t1,%lcoef2,%l_t2); lamda=2.; call print('Two test Cases - Penalized & Lasso formula to get Ridge':); call lasso2(%y,%x,%coef,%lcoef1,%l_t1,lamda,lresid1,1,0.); call lasso2(%y,%x,%coef,%lcoef2,%l_t2,lamda,lresid2,3,2.); call tabulate(%names,%lag,%coef,%se,%t,%lcoef1,%l_t1,%lcoef2,%l_t2); b34srun; %b34sif(&dorats.ne.0)then; b34sexec options open('rats.dat') unit(28) disp=unknown$ b34srun$ b34sexec options open('rats.in') unit(29) disp=unknown$ b34srun$ b34sexec options clean(28)$ b34srun$ b34sexec options clean(29)$ b34srun$ /; Uses logic from RATS User's Guide Version 6.1 Page 192 b34sexec pgmcall$ rats passasts pcomments('* ', '* Data passed from B34S(r) system to RATS', '* ', "display @1 %dateandtime() @33 ' Rats Version ' %ratsversion()" '* ') $ PGMCARDS$ * * Non centered results * cmoment # constant gasout{1 to 6} gasin{1 to 6} gasout linreg(cmoment) gasout # constant gasout{1 to 6} gasin{1 to 6} do row=1,13 compute %cmom(row,row)=%cmom(row,row)+2 end do row linreg(cmoment) gasout # constant gasout{1 to 6} gasin{1 to 6} b34sreturn$ b34srun $ b34sexec options close(28)$ b34srun$ b34sexec options close(29)$ b34srun$ b34sexec options /$ dodos(' rats386 rats.in rats.out ') dodos('start /w /r rats32s rats.in /run') dounix('rats rats.in rats.out')$ B34SRUN$ b34sexec options npageout WRITEOUT('Output from RATS',' ',' ') COPYFOUT('rats.out') dodos('ERASE rats.in','ERASE rats.out','ERASE rats.dat') dounix('rm rats.in','rm rats.out','rm rats.dat') $ B34SRUN$ %b34sendif; == ==LTS Least Trimmed Squares /; /; For further detail see Stokes (200x, Chapter 10) /; b34sexec options ginclude('gas.b34'); b34srun; b34sexec matrix; call loaddata; call load(lts :staging); call echooff; n=6; p=.8; iprint=1; call olsq(gasout gasin{1 to n} gasout{1 to n} :print :savex); call lts(%y,%x,%res,%names,%lag,%coef,%se,%t,p,newy,newx,ihold,iprint); b34srun; == ==LTS_1 Least Trimmed Squares OLS, MARS and GAM /; /; Application of LTS to MARS and GAM Models /; Is nonlinearity is Gas Data removed by outliers being removed? /; Or is detection made difficult? /; b34sexec options ginclude('gas.b34'); b34srun; b34sexec matrix; call loaddata; call load(lts :staging); call load(upmatrix :staging); call load(catname :staging); call echooff; n=6; p=.8; iprint=1; /; Before truncation what is nonlinearity? call olsq(gasout gasin{1 to n} gasout{1 to n} :print :savex); call marspline(gasout gasin{1 to n} gasout{1 to n} :print :mi 2 :nk 20); call gamfit(gasout gasin[predictor,3]{1 to n} gasout[predictor,3]{1 to n} :print); /; Now we truncate the data call olsq( gasout gasin{1 to n} gasout{1 to n} :savex); call lts(%y,%x,%res,%names,%lag,%coef,%se,%t,p,newy,newx,ihold, iprint); call marspline(newy newx :print :mi 2 :nk 20 :holdout ihold); call upmatrix(newx,%names ,%lag ,level(),newnames); add ='[predictor,3]'; call catname(newnames,add,augname,1); call gamfit(newy argument(augname) :holdout ihold :print); b34srun; == ==LTS_REC Recursive Least Trimmed Squares - Program Example /; /; Recursive Least Trimmed Squares /; b34sexec options ginclude('b34sdata.mac') member(res72); b34srun; b34sexec matrix; call loaddata; call load(lts :staging); call load(lts_rec :staging); call echooff; n_recur=12; p=.9; iprint=1; oldcoef=0; call olsq(lnq lnl lnk time :print :savex); call lts(%y,%x,%res,%names,%lag,%coef,%se,%t,p,newy,newx, ihold,iprint); call lts_rec; call print(%reccoef,%rect); b34srun; == ==LTS_REC2 Recursive Least Trimmed Squares - Subroutine Example b34sexec options ginclude('gas.b34'); b34srun; b34sexec matrix; call loaddata; call load(lts :staging); call load(lts_rec2 :staging); call load(lts_rec :staging); call echooff; /; Test running same command as a program if itestpgm=1 itestpgm=0; n=6; p=.8; iprint=1; call olsq(gasout gasin{1 to n} gasout{1 to n} :print :savex); call lts(%y,%x,%res,%names,%lag,%coef,%se,%t,p,newy,newx, ihold,iprint); oldcoef=1; n_recur=6; call lts_rec2(oldcoef,n_recur,p,newy,newx,%names,%lag,%coef,%se,%t, ihold,iprint); /; Test other way to run; if(itestpgm.ne.0)then; call print(' ':); call print('--------------------------------------':); call print(' Test of Program Version ':); call print('--------------------------------------':); call print(' ':); call olsq(gasout gasin{1 to n} gasout{1 to n} :print :savex); call lts(%y,%x,%res,%names,%lag,%coef,%se,%t,p,newy,newx, ihold,iprint); call lts_rec; endif; b34srun; == ==MARS_P Display a MARSPLINE Probit Model b34sexec options ginclude('b34sdata.mac') macro(murder)$ b34seend$ b34sexec matrix; call loaddata; call echooff; call load(mars_p :staging); yvar=afam(d1); /; this sets model call character(cc,'t y lf nw'); /; /; * Execute the initial Marspline model /; * Need to set yvar and cc /; * ----------------------------------- call marspline(yvar argument(cc) :probit :nk 20 :mi 2 :savebx :savemodel :df 2. :print); call mars_p; b34srun; == ==MCOVF Test Constant Covariance Matrix %b34slet dorats=0; b34sexec options ginclude('gas.b34'); b34srun; /; Illustrate mcovf function b34sexec matrix; call load(mcovf :staging); call loaddata; call echooff; call olsq(gasout gasin :savex); call print('test1a - usual case no lag',mcovf(%x,%res,0,0.0,0)); call print('test1b - usual case v ',mcovf(%x,%res,0,0.0,1)); call print('test2a - no residual ',mcovf(%x,0.0 ,0,0.0,0)); call print('test2b - no residual lag=3',mcovf(%x,0.0 ,3,0.0,0)); call print('test3 - lag = 1 ',mcovf(%x,%res,1,0.0,0)); call print('test4 - lag = 3 ',mcovf(%x,%res,3,0.0,0)); call print('test5 - lag = 2 damp=1. ',mcovf(%x,%res,2,1.0,0)); /; Optionally save workspace to Speakeasy /; case1a=mcovf(%x,%res,0,0.0,0); /; case1b=mcovf(%x,%res,0,0.0,1); /; case2a=mcovf(%x,0.0 ,0 0.0,0) /; case2b=mcovf(%x,0.0 ,3 0.0,0) /; case3=mcovf(%x,0.0 ,0,0.0,0); g /; case4=mcovf(%x,%res,1,0.0,0); /; case5=mcovf(%x,%res,2,1.0,0); /; call checkpoint; b34srun; %b34sif(&dorats.ne.0)%then; b34sexec options open('rats.dat') unit(28) disp=unknown$ b34srun$ b34sexec options open('rats.in') unit(29) disp=unknown$ b34srun$ b34sexec options clean(28)$ b34srun$ b34sexec options clean(29)$ b34srun$ b34sexec pgmcall$ rats passasts pcomments('* ', '* Data passed from B34S(r) system to RATS', '* ', "display @1 %dateandtime() @33 ' Rats Version ' %ratsversion()" '* ') $ PGMCARDS$ * linreg gasout / resids # gasin constant display 'Default case' mcov(damp=0.0,lags=0,print) * * resids # gasin constant display 'Nosquare Case' mcov(damp=1.0,print,nosquare) * * resids # constant gasin display 'Case without residuals lag=0' mcov(damp=0.0,lags=0,print) * * # constant gasin display 'Case without residuals but lag=3' mcov(damp=0.0,lags=3,print) * * # constant gasin display 'Case with residuals lag=1' mcov(damp=0.0,lags=1,print) * * resids # constant gasin display 'Case with residuals lag=3' mcov(damp=0.0,lags=3,print) * * resids # constant gasin display 'Case with residuals lag=2 Damp=1.' mcov(damp=1.0,lags=2,print) * * resids # constant gasin b34sreturn$ b34srun $ b34sexec options close(28)$ b34srun$ b34sexec options close(29)$ b34srun$ b34sexec options /$ dodos(' rats386 rats.in rats.out ') dodos('start /w /r rats32s rats.in /run') dounix('rats rats.in rats.out')$ B34SRUN$ b34sexec options npageout WRITEOUT('Output from RATS',' ',' ') COPYFOUT('rats.out') dodos('ERASE rats.in','ERASE rats.out','ERASE rats.dat') dounix('rm rats.in','rm rats.out','rm rats.dat') $ B34SRUN$ %b34sendif; == ==MARSDIAG Diagnostic Tests on a MARSPLINE Model b34sexec options ginclude('b34sdata.mac') member(gas); b34srun; b34sexec matrix; call loaddata; /; call echooff; call load(marsdiag :staging); call load(marsinfo :staging); m6=1; mm6=2; call olsq(gasout gasin{m6 to mm6} gasout{1 to mm6} :print); %res_ols = %res; %yhat_ols = %yhat; call marspline(gasout gasin{m6 to mm6} gasout{1 to mm6} :print :mi 2 :nk 20 :xx); %res_m2 = %res; %yhat_m2 = %yhat; call marsinfo; call marsdiag(%xx,c_sums,r_sums,1,1,' '); call marsdiag(%xx,c_sums,r_sums,2,2,'test.wmf'); call marsdiag(%xx,c_sums,r_sums,3,0,' '); /; call tabulate(%res_ols, %res_m2); /; call tabulate(%yhat_ols, %yhat_m2); call graph(%res_ols, %res_m2 :nolabel :nocontact :pgborder :heading 'MARS vs OLS residuals'); b34srun; == ==MARSVPLT Mars Variable Plot b34sexec options ginclude('b34sdata.mac') member(gas); b34srun; b34sexec matrix; call loaddata; /; call echooff; call load(marsdiag :staging); call load(marsvplt :staging); call load(marsinfo :staging); m6=1; mm6=2; call olsq(gasout gasin{m6 to mm6} gasout{1 to mm6} :print); %res_ols = %res; %yhat_ols = %yhat; call marspline(gasout gasin{m6 to mm6} gasout{1 to mm6} :print :mi 2 :nk 20 :xx); %res_m2 = %res; %yhat_m2 = %yhat; call marsdiag(%xx,c_sums,r_sums,1,1,' '); call marsdiag(%xx,c_sums,r_sums,2,2,'test.wmf'); call marsdiag(%xx,c_sums,r_sums,3,0,' '); call marsinfo; call marsvplt(%xx,'dd'); /; call tabulate(%res_ols, %res_m2); /; call tabulate(%yhat_ols, %yhat_m2); call graph(%res_ols,%res_m2 :nolabel :pgborder :nocontact); call print(%xx); b34srun; == ==MARSINFO GCV Information b34sexec options ginclude('b34sdata.mac') member(gas); b34srun; b34sexec matrix; call loaddata; call load(dispmars :staging); call load(marsINFO :staging); call echooff; call olsq(gasout gasin{1 to 6} gasout{1 to 6} :print); %res_ols = %res; %yhat_ols = %yhat; call mars(gasout gasin{1 to 6} gasout{1 to 6} :print :mi 2 :nk 20); %res_m = %res; %yhat_m = %yhat; call marspline(gasout gasin{1 to 6} gasout{1 to 6} :print :mi 6 :nk 35); %res_m2 = %res; %yhat_m2 = %yhat; call marsinfo; call graph(%res_ols, %res_m, %res_m2 :nolabel); b34srun; == ==NORM Matrix NORM & Condition /; Illustrates norm and cond /; Optionally Matlab can be called /; B34S tracks 100% Matlab / LAPACK /; %b34slet domatlab=0; b34sexec options ginclude('gas.b34'); b34srun; b34sexec matrix; call loaddata; call load(cov :staging); call load(cor :staging); call load(norm :staging); call load(cond :staging); call echooff; lag=3; /; call olsq(gasout gasout{1 to lag} gasin{1 to lag} :print :savex); call print('Tests on Real*8 Cholesky':); call print('++++++++++++++++++++++++':); xpx=transpose(%x)*%x; call print('Norm 1 xpx ',norm(xpx,1):); call print('Norm 2 xpx ',norm(xpx,2):); call print('Norm 3 xpx ',norm(xpx,3):); call print('Norm 4 xpx ',norm(xpx,4):); call print('Cond 1 xpx ',cond(xpx,1):); call print('Cond 2 xpx ',cond(xpx,2):); call print('Cond 3 xpx ',cond(xpx,3):); call print('Cond 4 xpx ',cond(xpx,4):); call print('Norm 1 x ',norm(%x,1):); call print('Norm 2 x ',norm(%x,2):); call print('Norm 3 x ',norm(%x,3):); call print('Norm 4 x ',norm(%x,4):); call print('Cond 2 x ',cond(%x,2):); test=inv(xpx,testcond :gmat); testcond=1.0/testcond; call print('Condition from LAPACK ',testcond:); test=inv(xpx,testcond ); testcond=1.0/testcond; call print('Condition from LINPACK ',testcond:); %b34sif(&domatlab.ne.0)%then; call makematlab(%x :file 'xdata.m'); call load(rmatlab); call rmatlab; pgmcards; format long g x=getb34s('xdata.m'); xpx=x'*x; % ' disp('norm(xpx,1)') disp( norm(xpx,1) ) disp('norm(xpx,2)') disp( norm(xpx,2) ) disp('norm(xpx,inf)') disp( norm(xpx,inf) ) disp('norm(xpx,"fro")') disp( norm(xpx,'fro')) disp('cond(xpx,1)') disp( cond(xpx,1) ) disp('cond(xpx,2)') disp( cond(xpx,2) ) disp('cond(xpx,inf)') disp( cond(xpx,inf) ) disp('cond(xpx,"fro")') disp( cond(xpx,'fro') ) disp('norm(x,1)') disp( norm(x,1) ) disp('norm(x,2)') disp( norm(x,2) ) disp('norm(x,inf)') disp( norm(x,inf) ) disp('norm(x,"fro")') disp( norm(x,'fro') ) quit b34sreturn; %b34sendif; b34srun; == ==NTOKIN Counts Number of tokins in a string b34sexec matrix; call character(cc,'10. 11 test y(10) jj=44 print'); call echooff; call load(ntokin :staging); call ntokin(cc,nfind,0,ibad); call print('All done tokenizing'); call print('# found was ',nfind); b34srun; == ==OLS_CON Constrained OLS Estimation b34sexec matrix; /; /; Constraint is of the form r=R*Bstar /; Bstar is the new coef vector. /; Bstar=Con_b defined below. /; We assume b1+b2+b3=1 and b1=-b2 /; /; Example where we use Theil (1971 page 44) to illustrate /; Constrained Least squares. OLS Coef should be 1.0 /; call load(ols_con :staging); n=10000; k=3; x=rn(matrix(n,k:)); beta=matrix(k,1:); beta(,1)=1.0; y=10.*rn(matrix(n,1:))+x*beta; /$ Constraint is where all coef sum to 1.0 /$ coef 2 = -1 * coef 1 call olsq(y x :noint :savex :print); small_r=matrix(2,1:1.0,0.0); r =matrix(2,3:1. 1. 1. 1.,1. 0.); call ols_con(%y,%x,%coef,con_b,small_r,r,ols_e,con_e); call print(con_b sumsq(ols_e),sumsq(con_e)); /; call tabulate(y,ols_e,con_e); b34srun; == ==POISSON Estimation of a Poisson Model /; Poisson test case /; routine poiss_SE SE1 & SE3 replicate LIMDEP & Rats for SE1 & SE3 %b34slet dolimdep=1; %b34slet dorats =1; %b34slet dob34s =1; b34sexec data heading('Example 2 Oster & Hilbe'); * Poisson model resp = f(x1,x2) ; * constant and SE should be -3.8103 and 1.7474 ; * See American Statistician Vol 62 Number 1 February 2008 "An Examination of Statistical Software Oackages for Parametric and Nonparametric Data Analysis Using Exact Methods" by Robert Oster and Joseph Hilbe; input resp x1 x2; datacards; 0 15 3 3 18 4 0 12 2 0 10 1 2 17 4 4 13 2 0 14 3 4 12 1 0 11 2 3 18 4 0 17 3 7 19 4 0 16 3 5 17 4 0 15 2 4 12 1 0 11 1 b34sreturn; b34srun; /$ b34sexec options open('limdep.in') unit(29) disp=unknown$ b34srun$ b34sexec options clean(29)$ b34srun$ %b34sif(&dob34s.eq.1)%then; b34sexec matrix; call loaddata; call load(poisson :staging); call load(poiss_se :staging); call print(poisson); call print(poiss_se); nseries=3; parm=vector(nseries:); rvec=vector(nseries:)+.01; y=vfam(resp); x=matrix(norows(y),nseries:); x(,1) =1.0; x(,2) =vfam(x1); x(,3) =vfam(x2); yfact=(-1.) * dlog(fact(resp)); syfact=sum(yfact); call print(syfact); call echooff; rvec=vector(nseries:)+.01; /; Starting from OLS call olsq(resp x1 x2 :print); rvec(1)=%coef(3); rvec(2)=%coef(1); rvec(3)=%coef(2); /; /; alternate start that is slower but will replicate RATS /; rvec=vector(nseries:)+.01; /$ /$ Note different SE's /$ call maxf2(func :name poisson :parms parm :maxit 3000 :ivalue rvec :print); call poiss_se(x,y,exp(xb),parm,se1,se2,se3,t1,t2,t3,1); /$ /$ start from same position /$ ll=vector(nseries:)-.1e+4; uu=vector(nseries:)+.1e+4; call cmaxf2(func :name poisson :parms parm :maxit 3000 :ivalue rvec :print :lower ll :upper uu); call poiss_se(x,y,exp(xb),parm,se1,se2,se3,t1,t2,t3,1); b34srun; %b34sendif; %b34sif(&dolimdep.eq.1)%then; b34sexec options noheader; b34srun; b34sexec pgmcall$ limdep ibat necho title('Results from limdep') pcomments('open ; output = limdep.out $') $ /$ pgmcards$ poisson ;lhs = resp ;rhs=one,x1,x2 $ b34sreturn$ b34srun $ b34sexec options close(29)$ b34srun$ b34sexec options position(1); b34srun; b34sexec options position(6); b34srun; b34sexec options dodos('limdep limdep.in') dounix('limdep < limdep.in'); b34srun b34sexec options npageout writeout('Output from LIMDEP',' ',' ') copyfout('limdep.out') writelog('Output from LIMDEP log') copyflog('trace.lim') dounix('rm limdep.in ', 'rm limdep.out', 'rm trace.lim ') dodos('erase limdep.in ', 'erase limdep.out', 'erase trace.lim ','cls') $ b34srun$ %b34sendif; %b34sif(&dorats.ne.0)%then; B34SEXEC OPTIONS OPEN('rats.dat') UNIT(28) DISP=UNKNOWN$ B34SRUN$ B34SEXEC OPTIONS OPEN('rats.in') UNIT(29) DISP=UNKNOWN$ B34SRUN$ B34SEXEC OPTIONS CLEAN(28)$ B34SRUN$ B34SEXEC OPTIONS CLEAN(29)$ B34SRUN$ B34SEXEC PGMCALL$ RATS PASSASTS pcomments('* ', '* Data passed from B34S(r) system to RATS', '* ', "display @1 %dateandtime() @33 ' Rats Version ' %ratsversion()" '* ') $ PGMCARDS$ * ddv(type=count) resp # constant x1 x2 ddv(type=count,robusterrors) resp # constant x1 x2 B34SRETURN$ B34SRUN $ B34SEXEC OPTIONS CLOSE(28)$ B34SRUN$ B34SEXEC OPTIONS CLOSE(29)$ B34SRUN$ B34SEXEC OPTIONS /$ dodos('start /w /r rats386 rats.in rats.out ') dodos('start /w /r rats32s rats.in /run') dounix('rats rats.in rats.out')$ B34SRUN$ B34SEXEC OPTIONS NPAGEOUT WRITEOUT('Output from RATS',' ',' ') COPYFOUT('rats.out') dodos('erase rats.in','erase rats.out','erase rats.dat') dounix('rm rats.in','rm rats.out','rm rats.dat') $ B34SRUN$ == ==POISS_SE Alternative Poisson SE's /; Poisson test case /; routine poiss_SE SE1 & SE3 replicate LIMDEP & Rats for SE1 & SE3 %b34slet dolimdep=1; %b34slet dorats =1; %b34slet dob34s =1; b34sexec data heading('Example 2 Oster & Hilbe'); * Poisson model resp = f(x1,x2) ; * constant and SE should be -3.8103 and 1.7474 ; * See American Statistician Vol 62 Number 1 February 2008 "An Examination of Statistical Software Oackages for Parametric and Nonparametric Data Analysis Using Exact Methods" by Robert Oster and Joseph Hilbe; input resp x1 x2; datacards; 0 15 3 3 18 4 0 12 2 0 10 1 2 17 4 4 13 2 0 14 3 4 12 1 0 11 2 3 18 4 0 17 3 7 19 4 0 16 3 5 17 4 0 15 2 4 12 1 0 11 1 b34sreturn; b34srun; /$ b34sexec options open('limdep.in') unit(29) disp=unknown$ b34srun$ b34sexec options clean(29)$ b34srun$ %b34sif(&dob34s.eq.1)%then; b34sexec matrix; call loaddata; call load(poisson :staging); call load(poiss_se :staging); call print(poisson); call print(poiss_se); nseries=3; parm=vector(nseries:); rvec=vector(nseries:)+.01; y=vfam(resp); x=matrix(norows(y),nseries:); x(,1) =1.0; x(,2) =vfam(x1); x(,3) =vfam(x2); yfact=(-1.) * dlog(fact(resp)); syfact=sum(yfact); call print(syfact); call echooff; rvec=vector(nseries:)+.01; /; Starting from OLS call olsq(resp x1 x2 :print); rvec(1)=%coef(3); rvec(2)=%coef(1); rvec(3)=%coef(2); /; /; alternate start that is slower but will replicate RATS /; rvec=vector(nseries:)+.01; /$ /$ Note different SE's /$ call maxf2(func :name poisson :parms parm :maxit 3000 :ivalue rvec :print); call poiss_se(x,y,exp(xb),parm,se1,se2,se3,t1,t2,t3,1); /$ /$ start from same position /$ ll=vector(nseries:)-.1e+4; uu=vector(nseries:)+.1e+4; call cmaxf2(func :name poisson :parms parm :maxit 3000 :ivalue rvec :print :lower ll :upper uu); call poiss_se(x,y,exp(xb),parm,se1,se2,se3,t1,t2,t3,1); b34srun; %b34sendif; %b34sif(&dolimdep.eq.1)%then; b34sexec options noheader; b34srun; b34sexec pgmcall$ limdep ibat necho title('Results from limdep') pcomments('open ; output = limdep.out $') $ /$ pgmcards$ poisson ;lhs = resp ;rhs=one,x1,x2 $ b34sreturn$ b34srun $ b34sexec options close(29)$ b34srun$ b34sexec options position(1); b34srun; b34sexec options position(6); b34srun; b34sexec options dodos('limdep limdep.in') dounix('limdep < limdep.in'); b34srun b34sexec options npageout writeout('Output from LIMDEP',' ',' ') copyfout('limdep.out') writelog('Output from LIMDEP log') copyflog('trace.lim') dounix('rm limdep.in ', 'rm limdep.out', 'rm trace.lim ') dodos('erase limdep.in ', 'erase limdep.out', 'erase trace.lim ','cls') $ b34srun$ %b34sendif; %b34sif(&dorats.ne.0)%then; B34SEXEC OPTIONS OPEN('rats.dat') UNIT(28) DISP=UNKNOWN$ B34SRUN$ B34SEXEC OPTIONS OPEN('rats.in') UNIT(29) DISP=UNKNOWN$ B34SRUN$ B34SEXEC OPTIONS CLEAN(28)$ B34SRUN$ B34SEXEC OPTIONS CLEAN(29)$ B34SRUN$ B34SEXEC PGMCALL$ RATS PASSASTS pcomments('* ', '* Data passed from B34S(r) system to RATS', '* ', "display @1 %dateandtime() @33 ' Rats Version ' %ratsversion()" '* ') $ PGMCARDS$ * ddv(type=count) resp # constant x1 x2 ddv(type=count,robusterrors) resp # constant x1 x2 B34SRETURN$ B34SRUN $ B34SEXEC OPTIONS CLOSE(28)$ B34SRUN$ B34SEXEC OPTIONS CLOSE(29)$ B34SRUN$ B34SEXEC OPTIONS /$ dodos('start /w /r rats386 rats.in rats.out ') dodos('start /w /r rats32s rats.in /run') dounix('rats rats.in rats.out')$ B34SRUN$ B34SEXEC OPTIONS NPAGEOUT WRITEOUT('Output from RATS',' ',' ') COPYFOUT('rats.out') dodos('erase rats.in','erase rats.out','erase rats.dat') dounix('rm rats.in','rm rats.out','rm rats.dat') $ B34S == ==P_L_EST Probit Logit estimation /; Test job /; dob34s => uses probit and loglin /; loglin => 2 times * -1. %b34slet dob34s =0$ %b34slet dorats =0$ %b34slet domatrix=1; b34sexec options ginclude('b34sdata.mac') macro(murder)$ b34seend$ %b34sif(&dob34s.ne.0)%then$ B34SEXEC PROBIT tola=.1e-14 $ MODEL D1 = T Y LF NW; B34Srun$ b34sexec loglin $ model d1 = T y lf nw; b34srun; B34SEXEC MLOGLIN IP=1 $ MODEL D1= T Y LF NW $ LEVEL D1(HAVELAW)$ B34SEEND$ %b34sendif$ %b34sif(&dorats.ne.0)%then$ b34sexec options ginclude('b34sdata.mac') macro(murder)$ b34seend$ b34sexec data set; rename t=tt; rename n=nn; b34srun; B34SEXEC OPTIONS HEADER$ B34SRUN$ b34sexec options open('rats.dat') unit(28) disp=unknown$ b34srun$ b34sexec options open('rats.in') unit(29) disp=unknown$ b34srun$ b34sexec options clean(28)$ b34srun$ b34sexec options clean(29)$ b34srun$ b34sexec pgmcall$ rats passasts PCOMMENTS('* ', '* Data passed from B34S(r) system to RATS', '* ') $ PGMCARDS$ * smpl lgt d1 # constant Tt Y LF NW prb d1 # constant Tt Y LF NW b34sreturn$ b34srun $ b34sexec options close(28)$ b34srun$ b34sexec options close(29)$ b34srun$ b34sexec options dodos('start /w /r rats32s rats.in /run') dounix('rats rats.in rats.out')$ B34SRUN$ b34sexec options npageout WRITEOUT('Output from RATS',' ',' ') COPYFOUT('rats.out') dodos('ERASE rats.in','ERASE rats.out','ERASE rats.dat') dounix('rm rats.in','rm rats.out','rm rats.dat') $ B34SRUN$ %b34sendif$ %b34sif(&domatrix.ne.0)%then; b34sexec options ginclude('b34sdata.mac') macro(murder)$ b34seend$ b34sexec matrix; call loaddata; call load(p_l_est :staging); call load(tlogit :staging); call echooff; %yin=d1; %xin=mfam(catcol(constant t y lf nw)); call print(%xin); call p_l_est('probit',%yin,%xin,%func,%coef,%se,%t,%yhat,'print',0); call print('Probit model':); call tabulate(%coef,%se,%t); call tabulate(%yin,%yhat); upper=.5; lower=.5; iprint=1; call character(cc,'Tests on probit Model'); call tlogit(%yin,%yhat,upper,lower,cc,ntruer,ntruep nfalser,nfalsep,nunclear,ptruer,pfalser,iprint); call p_l_est('logit', %yin,%xin,%func,%coef,%se,%t,%yhat,'print',0); call print('Logit model':); call tabulate(%coef,%se,%t); call tabulate(%yin,%yhat); upper=.5; lower=.5; iprint=1; call character(cc,'Tests on Logit Model'); call tlogit(%yin,%yhat,upper,lower,cc,ntruer,ntruep nfalser,nfalsep,nunclear,ptruer,pfalser,iprint); call print(' ':); call print(' Results with max2 ******************************':); call p_l_est('probit',%yin,%xin,%func,%coef,%se,%t,%yhat,'print',1); call print('Probit model':); call tabulate(%coef,%se,%t); call tabulate(%yin,%yhat); upper=.5; lower=.5; iprint=1; call character(cc,'Tests on probit Model'); call tlogit(%yin,%yhat,upper,lower,cc,ntruer,ntruep nfalser,nfalsep,nunclear,ptruer,pfalser,iprint); call p_l_est('logit', %yin,%xin,%func,%coef,%se,%t,%yhat,'print',1); call print('Logit model':); call tabulate(%coef,%se,%t); call tabulate(%yin,%yhat); upper=.5; lower=.5; iprint=1; call character(cc,'Tests on Logit Model'); call tlogit(%yin,%yhat,upper,lower,cc,ntruer,ntruep nfalser,nfalsep,nunclear,ptruer,pfalser,iprint); b34srun; %b34sendif; == ==PANEL_LIB Test of Panel Subroutine Library /; /; Panel subroutines tests. /; panel_t => transpose panel /; panel_df => Difference panel /; panel_fe => Removes Mean used with one-way fixed effects /; panel2fe => Removes Means used with two-way fixed effects /; can do time and individual effects without transpose. /; pfe_1way => Driver routine for one way Fixed Effects /; pfe_2way => Driver routine for two way Fixed effects /; /; Uses RATS to validate one and two way tests /; %b34slet runrats=1; b34sexec options ginclude('panel_data.mac') member(grunfeld); b34srun; /; b34sexec list; b34srun; /$ Shows ECOMP and REG %b34sif(&runrats.ne.0)%then; b34sexec options open('rats.dat') unit(28) disp=unknown$ b34srun$ b34sexec options open('rats.in') unit(29) disp=unknown$ b34srun$ b34sexec options clean(28)$ b34srun$ b34sexec options clean(29)$ b34srun$ b34sexec pgmcall$ rats passasts pcomments('* ', '* Data passed from B34S(r) system to RATS', '* ', "display @1 %dateandtime() @33 ' Rats Version ' %ratsversion()" '* ') $ PGMCARDS$ * cal(panel=20) 1935 1 1 * print linreg invest / resids # constant f c pregress(method=fixed,effects=both) invest # f c pregress(method=fixed,effects=individual) invest # f c pregress(method=fixed,effects=time) invest # f c pregress(method=fixed,effects=both) invest # constant f c pregress(method=random,effects=individual) invest # constant f c pregress(method=random,effects=time) invest # constant f c pregress(method=fd,effects=individual) invest # constant f c pregress(method=fd,effects=time) invest # constant f c * pregress(method=sur) invest * # constant f c b34sreturn$ b34srun $ b34sexec options close(28)$ b34srun$ b34sexec options close(29)$ b34srun$ b34sexec options /$ dodos(' rats386 rats.in rats.out ') dodos('start /w /r rats32s rats.in /run') dounix('rats rats.in rats.out')$ B34SRUN$ b34sexec options npageout WRITEOUT('Output from RATS',' ',' ') COPYFOUT('rats.out') dodos('ERASE rats.in','ERASE rats.out','ERASE rats.dat') dounix('rm rats.in','rm rats.out','rm rats.dat') $ B34SRUN$ %b34sendif; b34sexec reg $ model invest = f c$ b34srun; b34sexec ecomp regfirst bothp nreg=10 nper=20$ model invest = f c$ b34srun$ b34sexec ecomp regfirst iswith bothp nreg=10 nper=20$ model invest = f c$ b34srun$ b34sexec matrix; call loaddata; call load(panel_lib :staging); call olsq(invest f c :print :savex); %h_names=%names(integers(1,%k-1)); /; /; hold data /; call echooff; %xx=%x; %yy=%y; /; call names; itest1=1; itest2=1; itest3=1; itest4=1; if(itest1.eq.1)then; call print(' ':); call print('+++++++++++++++++++++++++++++++++++++++++++++++':); call print(' ':); call panel2fe(%yy,%xx,10,20,%ynew,%xnew,2); call print('Two Way Fixed Effects Estimator panel2fe=2':); call deletecol(%xnew,nocols(%xnew)); call olsq(%ynew %xnew :print :noint); endif; if(itest2.eq.1)then; call print(' ':); call print('+++++++++++++++++++++++++++++++++++++++++++++++':); call print(' ':); call print('Data sorted Before Time fixed effect':); call panel_t( %yy,%xx,10,20,%ynew1,%xnew1); call panel2fe(%ynew1,%xnew1,20,10,%ynew,%xnew,0); call print('Fixed Effects Estimator time correction panel2fe=0':); call deletecol(%xnew,nocols(%xnew)); call olsq(%ynew %xnew :print :noint); call print(' ':); call print('+++++++++++++++++++++++++++++++++++++++++++++++':); call print(' ':); call print('Data not sorted Time Fixed Effect. panel_t not used':); call panel2fe(%yy,%xx,10,20,%ynew,%xnew,1); call print('Fixed Effects Estimator time correction panel2fe=1':); call deletecol(%xnew,nocols(%xnew)); call olsq(%ynew %xnew :print :noint); endif; if(itest3.eq.1)then; call print(' ':); call print('+++++++++++++++++++++++++++++++++++++++++++++++':); call print(' ':); call print('Using panel2fe code':); call panel2fe(%yy,%xx,10,20,%ynew,%xnew,0); call print('Fixed Effects Estimator individual lsvd':); call deletecol(%xnew,nocols(%xnew)); call olsq(%ynew %xnew :print :noint); call print(' ':); call print('+++++++++++++++++++++++++++++++++++++++++++++++':); call print(' ':); call print('Using panel_fe code':); call panel_fe(%yy,%xx,10,20,%ynew,%xnew); call print('Fixed Effects Estimator individual lsvd':); call deletecol(%xnew,nocols(%xnew)); call olsq(%ynew %xnew :print :noint); endif; if(itest4.ne.0)then; call print(' ':); call print('+++++++++++++++++++++++++++++++++++++++++++++++':); call print(' ':); call print('Differencing of a n by t panel':); call panel_df(%yy,%xx,10,20,%ynew,%xnew); call print('FD Fixed Effects Estimator':); call deletecol(%xnew,nocols(%xnew)); call olsq(%ynew %xnew :print); call olsq(%ynew %xnew :print :noint); endif; b34srun; == ==PANEL_LIB2 Test of Panel Subroutine Library /; /; Panel subroutines /; panel_t => transpose panel /; panel_df => Diffference panel /; panel_fe => Removes Mean /; panel2fe => Two Way panel Model /; b34sexec options ginclude('panel_data.mac') member(grunfeld); b34srun; /; b34sexec list; b34srun; /$ Shows ECOMP and REG b34sexec reg $ model invest = f c$ b34srun; b34sexec ecomp regfirst bothp nreg=10 nper=20$ model invest = f c$ b34srun$ b34sexec ecomp regfirst iswith bothp nreg=10 nper=20$ model invest = f c$ b34srun$ b34sexec matrix; call loaddata; call load(panel_lib :staging); call olsq(invest f c :print :savex); /; /; hold data /; %xx=%x; %yy=%y; /; call names; call echooff; program test_p_t; /; /; Program to illustrate what is done with panel_t /; n=20; x =grid(1.,dfloat(n),1.); call print(x); y =rn(x); yy=rn(x); i=integers(1,5); do j=1,(n/5); yy(i)=dfloat(j); i=i+5; enddo; bigx=catcol(x,yy); call print('RAW data file that is 5,4':); call print(catcol(y,bigx)); call panel_t(y,bigx,4,5,ynew,xnew); call print(catcol(ynew,xnew)); call panel_t(ynew,xnew,5,4,yold,xold); call print(catcol(yold,xold)); return; end; /; call test_p_t; call print(' ':); call print('+++++++++++++++++++++++++++++++++++++++++++++++':); call print(' ':); call print('Transpose n by t panel':); call tabulate(fn,yr,invest,f,c); call panel_t( %yy,%xx,10,20,%ynew1,%xnew1); call tabulate(%yy,%ynew1); call print(%xx,%xnew1); call print(' ':); call print('+++++++++++++++++++++++++++++++++++++++++++++++':); call print(' ':); call print('Differencing of a n by t panel':); call panel_df(%yy,%xx,10,20,%ynew,%xnew); call print('FD Fixed Effects Estimator':); call deletecol(%xnew,nocols(%xnew)); call olsq(%ynew %xnew :print); call olsq(%ynew %xnew :print :noint); call print(' ':); call print('+++++++++++++++++++++++++++++++++++++++++++++++':); call print(' ':); call print('Data sorted other way':); call panel_t( %yy,%xx,10,20,%ynew1,%xnew1); call panel_df(%ynew1,%xnew1,20,10,%ynew,%xnew); call print('FD Fixed Effects Estimator':); call deletecol(%xnew,nocols(%xnew)); call olsq(%ynew %xnew :print); call olsq(%ynew %xnew :print :noint); call print(' ':); call print('+++++++++++++++++++++++++++++++++++++++++++++++':); call print(' ':); call print('mean removed of a n by t panel':); call panel_fe(%yy,%xx,10,20,%ynew,%xnew); call print('Fixed Effects Estimator':); call deletecol(%xnew,nocols(%xnew)); call olsq(%ynew %xnew :print :noint); call print(' ':); call print('+++++++++++++++++++++++++++++++++++++++++++++++':); call print(' ':); call print('Data sorted other way':); call panel_t( %yy,%xx,10,20,%ynew1,%xnew1); call panel_fe(%ynew1,%xnew1,20,10,%ynew,%xnew); call print('Fixed Effects Estimator':); call deletecol(%xnew,nocols(%xnew)); call olsq(%ynew %xnew :print :noint); b34srun; == ==PFE_1WAY Tests both one and two way Fixed Effects Models /; /; Panel subroutines in library. /; panel_t => transpose panel /; panel_df => Difference panel /; panel_fe => Removes Mean used with one-way fixed effects /; panel2fe => Removes Means used with two-way fixed effects /; can do time and individual effects without transpose. /; pfe_1way => Driver routine for one way Fixed Effects /; pfe_2way => Driver routine for two way Fixed effects /; /; Uses RATS to validate one and two way tests /; %b34slet runrats=0; b34sexec options ginclude('panel_data.mac') member(grunfeld); b34srun; /; b34sexec list; b34srun; b34sexec matrix; call echooff; call loaddata; call load(panel_lib :staging); call olsq(invest f c :print :savex); /; /; do individual, time/region and difference Fixed Effects /; /; n= # of individuals /; t = # of tiem periods /; n=10; t=20; call pfe_1way(%y,%x,%names,n,t,0); call pfe_1way(%y,%x,%names,n,t,1); call pfe_1way(%y,%x,%names,n,t,2); /; Do two way fixed effects call pfe_2way(%y,%x,%names,n,t); b34srun; /; This code calls RATS to validate calculations %b34sif(&runrats.ne.0)%then; b34sexec options open('rats.dat') unit(28) disp=unknown$ b34srun$ b34sexec options open('rats.in') unit(29) disp=unknown$ b34srun$ b34sexec options clean(28)$ b34srun$ b34sexec options clean(29)$ b34srun$ b34sexec pgmcall$ rats passasts pcomments('* ', '* Data passed from B34S(r) system to RATS', '* ', "display @1 %dateandtime() @33 ' Rats Version ' %ratsversion()" '* ') $ PGMCARDS$ * cal(panel=20) 1935 1 1 * print linreg invest / resids # constant f c pregress(method=fixed,effects=both) invest # f c pregress(method=fixed,effects=individual) invest # f c pregress(method=fixed,effects=time) invest # f c pregress(method=fixed,effects=both) invest # constant f c pregress(method=random,effects=individual) invest # constant f c pregress(method=random,effects=time) invest # constant f c pregress(method=random,effects=both) invest # constant f c pregress(method=fd,effects=individual) invest # constant f c pregress(method=fd,effects=time) invest # constant f c b34sreturn$ b34srun $ b34sexec options close(28)$ b34srun$ b34sexec options close(29)$ b34srun$ b34sexec options /$ dodos(' rats386 rats.in rats.out ') dodos('start /w /r rats32s rats.in /run') dounix('rats rats.in rats.out')$ B34SRUN$ b34sexec options npageout WRITEOUT('Output from RATS',' ',' ') COPYFOUT('rats.out') dodos('ERASE rats.in','ERASE rats.out','ERASE rats.dat') dounix('rm rats.in','rm rats.out','rm rats.dat') $ B34SRUN$ %b34sendif; == ==QR_SMALL Breaks appart QR Factorization b34sexec matrix; call load(qr_small :staging); /; Assume r=qr(x,q:); /; q1 =submatrix(q,1,norows(q),1,nocols(x)); /; q2 =submatrix(q,1,norows(q),nocols(x)+1, nocols(q)); /; rsmall=submatrix(r,1,nocols(x),1,nocols(x)); x=rn(matrix(10,3:)); call print('x(10,3) First r(3,3). Full r(10,3)':); call print(qr(x)); call print(qr(x:)); r1=qr(x,q1); call print('Economy R Q m>n => r(n,n) q(m,n)':); call print(r1,q1,x); call print('test of factorization ',transpose(r1)*r1,transpose(x)*x); call print('test of transpose(q)*q',transpose(q1)*q1); call print(' ':); call print('Full R Q m>n => r(m,n) q(m,m)':); r1=qr(x,q:); call print(r1,q,x); call print('test of factorization',transpose(r1)*r1,transpose(x)*x); call print('test of transpose(q)*q',transpose(q)*q); call print(' ':); call qr_small(x,q,r1,q1,q2,rsmall); call print(q,r1,q1,q2,rsmall); call print('Application to OLS'); nn=10; x=rn(matrix(nn,4:)); x(,1)=1.0; y=rn(vector(nn:)); y=to_vector(sumrows(x)+y); call olsq(y,x :noint :print); r=qr(x,q:); call qr_small(x,q,r,q1,q2,rsmall); call print('Test regression has coef = 1':); call print(r,q); yhat =to_vector(q1*transpose(q1)*y); error=to_vector(q2*transpose(q2)*y); call tabulate(y,%yhat,%res,yhat,error); b34srun; == ==QR_SMALL2 Error Buildup /$ Illustrates OLS Capability under Matrix Command b34sexec options ginclude('b34sdata.mac') member(gas); b34srun; b34sexec matrix; call load(qr_small :staging); call echooff; call loaddata; nlag=4; call olsq(gasout gasin{0 to nlag} gasout{1 to nlag} :print :savex); call print(ccf(%y, %res)); call print(ccf(%yhat,%res)); /; /; Large QR used for illustration. Q2 is large!!! /; error2 equation uses q1 only => uses the economy qr. This is the /; best way to proceed /; r=qr(%x,q:); call qr_small(%x,q,r,q1,q2,r_small); /; call print(q,q1,q2); yhat =q1*transpose(q1)*%y; error =q2*transpose(q2)*%y; error2=%y - yhat; beta=inv(r_small)*transpose(q1)*%y; call print('Beta from QR ',beta); call print(ccf(%y,error)); call print(ccf(yhat,error)); /; call tabulate(%y,%yhat,yhat,%res,error,error2); call print(' ':); call print('Study Error Buildup using Cholesky':); call print(' ':); /; excessive problem maxlag=40; chol_r=vector(maxlag:); qr_r =vector(maxlag:); do i=1,maxlag; /; /; :qr call uses linpack to get OLS. This is close to LAPACK QR( ) /; /; call olsq(gasout gasin{0 to i} gasout{1 to i} :qr :savex); call olsq(gasout gasin{0 to i} gasout{1 to i} :savex); /; /; Use economy size qr to save space!! /; r=qr(%x,q); qr_yhat =q*transpose(q)*%y; qr_error =%y - qr_yhat; chol_r(i)=ccf(%yhat , %res); qr_r(i) =ccf(qr_yhat,qr_error); enddo; call tabulate(chol_r,qr_r :title 'As maxlag increases accuracy declines'); b34srun; == ==RIDGE Ridge Regressioin /; /; For further detail see Stokes (200x, Chapter 10) /; %b34slet dorats=1; b34sexec options ginclude('gas.b34'); b34srun; b34sexec matrix; call loaddata; call load(ridge :staging); call echooff; k=6; call olsq(gasout gasout{1 to k} gasin{1 to k} :print :savex); lamda=2.0; call ridge(%y,%x,lamda,%coef,%names,%lag,ridge_c,0); call ridge(%y,%x,lamda,%coef,%names,%lag,ridge_c,1); b34srun; %b34sif(&dorats.ne.0)then; b34sexec options open('rats.dat') unit(28) disp=unknown$ b34srun$ b34sexec options open('rats.in') unit(29) disp=unknown$ b34srun$ b34sexec options clean(28)$ b34srun$ b34sexec options clean(29)$ b34srun$ /; Uses logic from RATS User's Guide Version 6.1 Page 192 b34sexec pgmcall$ rats passasts pcomments('* ', '* Data passed from B34S(r) system to RATS', '* ', "display @1 %dateandtime() @33 ' Rats Version ' %ratsversion()" '* ') $ PGMCARDS$ * * Non centered results * cmoment # constant gasout{1 to 6} gasin{1 to 6} gasout linreg(cmoment) gasout # constant gasout{1 to 6} gasin{1 to 6} do row=1,13 compute %cmom(row,row)=%cmom(row,row)+2 end do row linreg(cmoment) gasout # constant gasout{1 to 6} gasin{1 to 6} b34sreturn$ b34srun $ b34sexec options close(28)$ b34srun$ b34sexec options close(29)$ b34srun$ b34sexec options /$ dodos(' rats386 rats.in rats.out ') dodos('start /w /r rats32s rats.in /run') dounix('rats rats.in rats.out')$ B34SRUN$ b34sexec options npageout WRITEOUT('Output from RATS',' ',' ') COPYFOUT('rats.out') dodos('ERASE rats.in','ERASE rats.out','ERASE rats.dat') dounix('rm rats.in','rm rats.out','rm rats.dat') $ B34SRUN$ %b34sendif; == ==RNW_SE Newey West Standard Error /$ /$ Test Problem to replicate Newey-West Test output in Greene Edition 4 /$ %b34slet dorats=1; b34sexec options ginclude('greene.mac') member(a13_1); b34srun; b34sexec matrix; call loaddata; call load(rnw_se :staging); call load(mcovf :staging); call echooff; /$ Simple call call olsq(realnvst realgnp realint :white :savex :print); call echooff; call rnw_se; lag=4; damp=1.0; call echooff; /$ Direct call call nw_se(%names,%lag,%coef,%xpxinv,%res, damp,%x,%se,%whitese,%nwse,lag,nw,white,1); /$ See how sensitive to a range of lags!!!! do i=4,8; lag=i; damp=1.0; call print('lag was ',i); call nw_se(%names,%lag,%coef,%xpxinv,%res, damp,%x,%se,%whitese,%nwse,lag,nw,white,1); enddo; realnvst=r8tor16(realnvst); realgnp =r8tor16(realgnp); realint =r8tor16(realint); /$ Real*16 tests call olsq(realnvst realgnp realint :white :savex :print); call rnw_se; b34srun; %b34sif(&dorats.eq.1)%then; b34sexec options open('rats.dat') unit(28) disp=unknown$ b34srun$ b34sexec options open('rats.in') unit(29) disp=unknown$ b34srun$ b34sexec options clean(28)$ b34srun$ b34sexec options clean(29)$ b34srun$ b34sexec pgmcall$ rats passasts PCOMMENTS('* ', '* Data passed from B34S(r) system to RATS', '* ') $ PGMCARDS$ * linreg(robusterrors,damp=1.0,lag=4) realnvst * * # constant realgnp realint * linreg(robusterrors,damp=1.0,lag=8) realnvst * * # constant realgnp realint * White Statistic linreg(robusterrors) realnvst * * # constant realgnp realint b34sreturn$ b34srun $ b34sexec options close(28)$ b34srun$ b34sexec options close(29)$ b34srun$ b34sexec options /$ dodos(' rats386 rats.in rats.out ') dodos('start /w /r rats32s rats.in /run') dounix('rats rats.in rats.out')$ B34SRUN$ b34sexec options npageout WRITEOUT('Output from RATS',' ',' ') COPYFOUT('rats.out') dodos('ERASE rats.in','ERASE rats.out','ERASE rats.dat') dounix('rm rats.in','rm rats.out','rm rats.dat') $ B34SRUN$ %b34sendif; == ==SVD_OLS OLS Using SVD Approach b34sexec options ginclude('gas.b34'); b34srun; b34sexec matrix; call loaddata; call load(svd_ols :staging); call load(svd2_ols :staging); call print(svd_ols); call olsq(gasout gasin{1 to 6} gasout{1 to 6} :print :savex); call echooff; call svd_ols(%y,%x,u,v,s,beta_svd,se_svd,pc_coef,ss_svd,rss,ibad); call print(ss_svd); call names; call tabulate(%names,%lag,pc_coef,beta_svd,se_svd,%coef,%se); call svd2_ols(%y,%x,u,v,s,beta_svd,se_svd,pc_coef,ss_svd,rss,ibad); call print(ss_svd); call names; call tabulate(%names,%lag,pc_coef,beta_svd,se_svd,%coef,%se); call print('Real*16 Results ************************':); gasin16 =r8tor16(gasin); gasout16=r8tor16(gasout); call olsq(gasout16 gasin16{1 to 6} gasout16{1 to 6} :print :savex); call echooff; call svd_ols(%y,%x,u1,v,s,beta_svd,se_svd,pc_coef,ss_svd,rss,ibad); call print(ss_svd); call tabulate(%names,%lag,pc_coef, beta_svd,se_svd,%coef,%se); b34srun; == ==SVD2_OLS OLS Using SVD Approach b34sexec options ginclude('gas.b34'); b34srun; b34sexec matrix; call loaddata; call load(svd_ols :staging); call load(svd2_ols :staging); call print(svd_ols); call olsq(gasout gasin{1 to 6} gasout{1 to 6} :print :savex); call echooff; call svd_ols(%y,%x,u,v,s,beta_svd,se_svd,pc_coef,ss_svd,rss,ibad); call print(ss_svd); call names; call tabulate(%names,%lag,pc_coef,beta_svd,se_svd,%coef,%se); call svd2_ols(%y,%x,u,v,s,beta2svd,se2svd,pc2coef,ss2svd,rss,ibad); call print(ss_svd); call names; call tabulate(%names,%lag,pc2coef,beta2svd,se2svd,%coef,%se); b34srun; == ==TESTACC Test Accuracy using Residual Orthognality /; /; Charles Renfro Test Regression. /; Uses Real*4, Real*8, Real*16 and VPA to test the correlation /; Between the residual and the right hand sides of the equation. /; This is a simple case. The next problem is harder /; /; This test problem implements some of the ideas in the paper /;"Investigating the effect of dsata precisioon and OLS calculation /; methods on obtaining residuals that are orthogonal to the right /; hand side variables" by Houston H. Stokes %b34slet runmat1=0; %b34slet runmat2=0; %b34slet dovpa =0; b34sexec options ginclude('b34sdata.mac') member(grunfeld_4); b34srun; b34sexec matrix; call loaddata; call echooff; call load(cov :staging); call load(cor :staging); call load(norm :staging); call load(cond :staging); call load(testacc :staging); /; /; Set parameters for tests /; /; call print(' ':); call print('GE Equation':); call print('-----------':); call print(' ':); call print(' ':); call echooff; /; itest=2; call olsq(gei gef gec :print :savex); call testacc(%x,%res,itest); %yhold=%y; %xhold=%x; r8_coef=%coef; call olsq(gei gef gec :print :savex :qr); call testacc(%x,%res,itest); %ytest=r8tor16(%yhold); %xtest=r8tor16(%xhold); call olsq(%ytest %xtest :print :noint :savex); call testacc(%xtest,%res,itest); r16_coef=%coef; call olsq(%ytest %xtest :print :noint :savex :qr); call testacc(%xtest,%res,itest); call print('Real*4 Tests ':); call print('____________ ':); %ytest=r8tor4(%yhold); %xtest=r8tor4(%xhold); beta = inv(transpose(%xtest)*%xtest)*transpose(%xtest)*%ytest; call print('Beta from Real*4 using inverse ':); r4_coef=beta; call sigd(16); call print(r8_coef,r16_coef,r4_coef); call print(beta); i=nocols(%xtest); %res=%xtest*beta-%ytest; call testacc(%xtest,%res,itest); /; QR approach r4_r=qr( %xtest ,r4_q); /; r8_r=qr(r4tor8(%xtest),r8_q); /; call print('R looked at from real 4 to real 8 and from real 4', /; r8_r,r4_r); r4_yhat=r4_q*transpose(r4_q)*%ytest; r4_res =%ytest-r4_yhat; /; /; Beta not needed to get yhat from QR !! /; r4_beta=inv(r4_r)*transpose(r4_q)*%ytest; call print('Beta from Real*4 using QR ':); call print(r4_beta); call testacc(%xtest,r4_res,itest); /; Alternative less accurate way to get residual r4_res_a=%ytest-%xtest*r4_beta; call print('Correlations from Real*4 QR via Beta':); call testacc(%xtest,r4_res_a,itest); %b34sif(&dovpa.ne.0)%then; call print('VPA Tests ':); call print('_________ ':); %ytest=vpa(%yhold); %xtest=vpa(%xhold); beta = inv(transpose(%xtest)*%xtest)*transpose(%xtest)*%ytest; call print('beta from VPA ':); call print(beta); %res=%xtest*beta-%ytest; call testacc(%xtest,%res,itest); %b34sendif; %b34sif(&runmat1.ne.0)%then; call makematlab(%xhold :file 'xdata1.m'); call makematlab(%yhold :file 'ydata1.m'); call load(rmatlab); call rmatlab; call dodos('pause'); call dodos('copy test.m test1.m'); /$ /$ Matlab Commands /$ pgmcards; x=getb34s('xdata1.m'); y=getb34s('ydata1.m'); beta8=inv(x'*x)*x'*y; res8=y-x*beta8; yhat_8=x*beta8; newx8=x; newx8(:,3)=res8; disp('Matlab LU Real*8 results') c8=corr(newx8) beta4=inv(single(x)'*single(x))*single(x)'*single(y); yhat_4=single(x)*beta4; res4 =single(y)-single(x)*beta4; newx4=single(x); newx4(:,3)=res4; disp('Matlab LU real*4 results') c4=corr(newx4) betatest=[beta8 double(beta4) ] % % qr cond_x=cond(x); [q8,r8]=qr(x,0); [q4,r4]=qr(single(x),0); yhat_q8=q8*q8'*y; yhat_q4=q4*q4'*single(y); % plot(yhat_q8-yhat_q4) % [yhat_q4 yhat_q8 yhat_8 yhat_4] res_q8= y -yhat_q8; res_q4=single(y)-yhat_q4; % Get residual another way via beta to test beta_q4=inv(r4)*q4'*single(y); % ' yhat_q4_alt=single(x)*beta_q4; res_q4_alt=single(y)-yhat_q4_alt; newx8_q=x; newx4_q=single(x); newx8_q(:,3)=res_q8; newx4_q(:,3)=res_q4; c8_q=corr(newx8_q); c4_q=corr(newx4_q); disp('alt c4_q via beta'); newx4_q_alt=single(x); newx4_q_alt(:,3)=res_q4_alt; c4_q_alt=corr(newx4_q_alt); disp('Matlab QR real*8 and QR real*4') testres=[res8 res_q8 res4 res_q4 res_q4_alt ]; testaccr =[c8(3,:)' c8_q(3,:)' c4(3,:)' c4_q(3,:)' c4_q_alt(3,:)']'; testaccr % Next two commands pass Matlab data back to Matrix % makeb34s('e3.dat',e3); % makeb34s('evec3.dat',evec3); % quit % testres=[res8 res_q8 res4 res_q4] % Next two commands pass Matlab data back to Matrix % makeb34s('e3.dat',e3); % makeb34s('evec3.dat',evec3); disp('+++++++++++++ Matlab Ending ++++++++++++++++++') quit b34sreturn; %b34sendif; b34srun; /; b34sexec options ginclude('gas.b34'); b34srun; b34sexec matrix; call loaddata; call load(cov :staging); call load(cor :staging); call load(norm :staging); call load(cond :staging); call load(testacc :staging); call echooff; lag=3; /; /; This options takes space /; itest=2; dovpa=0; %b34sif(&dovpa.ne.0)%then; dovpa=1; %b34sendif; doreal4=1; call olsq(gasout gasout{1 to lag} gasin{1 to lag} :print :savex :diag); call print(%coef,%t); call print('Tests on Real*8 Cholesky':); call print('++++++++++++++++++++++++':); x_r8=%x; y_r8=%y; call testacc(x_r8,%res,0); call testacc(x_r8,%res,itest); call olsq(gasout gasout{1 to lag} gasin{1 to lag} :print :savex :qr :diag); call print('Tests on Real*8 QR':); call print('++++++++++++++++++':); call testacc(x_r8,%res,itest); %yhold=%y; %xhold=%x; xpx=transpose(%x)*%x; %ytest=r8tor16(%yhold); %xtest=r8tor16(%xhold); call olsq(%ytest %xtest :print :diag :noint :savex); call print('Tests on Real*16 Cholesky':); call print('+++++++++++++++++++++++++':); call testacc(%xtest,%res,itest); call olsq(%ytest %xtest :print :diag :noint :qr :savex); call print('Tests on Real*16 QR':); call testacc(%xtest,%res,itest); if(dovpa.ne.0)then; call print('VPA Tests ':); call print('_________ ':); %ytest=vpa(%yhold); %xtest=vpa(%xhold); beta = inv(transpose(%xtest)*%xtest)*transpose(%xtest)*%ytest; call print('beta from VPA ':); call print(beta); i=nocols(%xtest); %res=%xtest*beta-%ytest; call testacc(%xtest,%res,itest); endif; /; This example show how key method is when data is not /; saved with enough precision if(doreal4.ne.0)then; call print('Real*4 Tests ':); call print('_________ ':); %ytest=r8tor4(%yhold); %xtest=r8tor4(%xhold); beta = inv(transpose(%xtest)*%xtest)*transpose(%xtest)*%ytest; r4_res=%xtest*beta-%ytest; call print('beta from Real*4 using inverse ':); call print(beta); call testacc(%xtest,r4_res,itest); /; QR approach r4_r=qr( %xtest ,r4_q); r4_yhat=r4_q*transpose(r4_q)*%ytest; r4_res_q =%ytest-r4_yhat; /; /; Beta not needed to get yhat from QR !! /; r4_beta=inv(r4_r)*transpose(r4_q)*%ytest; call print('beta from Real*4 using QR ':); call print(r4_beta); call testacc(%xtest,r4_res_q,itest); /; Alternative less accurate way to get residual r4_res_a=%ytest-%xtest*r4_beta; call print('Correlations from Real*4 QR via Beta':); call testacc(%xtest,r4_res_a,itest); /; ############################################################### %b34sif(&runmat2.ne.0)%then; call makematlab(x_r8 :file 'xdata2.m'); call makematlab(y_r8 :file 'ydata2.m'); call load(rmatlab); call rmatlab; call dodos('copy test.m test2.m'); /$ Optional * call getmatlab(e3 :file 'e3.dat'); * call getmatlab(evec3 :file 'evec3.dat'); * call print('Matlab Answers',e3,evec3); /$ Running Matlab script under B34S Matrix /$ Tasks: 1. Pass Data from B34S to Matlab /$ 2. Do work in Matlab /$ 3. Bring data back from Matlab to B34S /$ not implemented /$ user needs to type quit in matlab to get /$ back to b34s /$ /$ Matlab Commands /$ pgmcards; x=getb34s('xdata2.m'); y=getb34s('ydata2.m'); beta8=inv(x'*x)*x'*y; res8=y-x*beta8; yhat_8=x*beta8; newx8=x; newx8(:,7)=res8; disp('Matlab LU real*8 results') c8=corr(newx8) beta4=inv(single(x)'*single(x))*single(x)'*single(y); disp('++++++++++++++++ beta4-beta8 +++++++++++++++++++++++') beta4-beta8 yhat_4=single(x)*beta4; res4 =single(y)-single(x)*beta4; newx4=single(x); newx4(:,7)=res4; disp('Matlab LU real*4 results') c4=corr(newx4) betatest=[beta8 double(beta4) ] % qr disp('Matlab QR real*8 and QR real*4') [q8,r8]=qr(x,0); [q4,r4]=qr(single(x),0); yhat_q8=q8*q8'*y; yhat_q4=q4*q4'*single(y); % plot(yhat_q8-yhat_q4) % [yhat_q4 yhat_q8 yhat_8 yhat_4] res_q8= y -yhat_q8; res_q4=single(y)-yhat_q4; % Get residual another way via beta to test beta_q4=inv(r4)*q4'*single(y); % ' yhat_q4_alt=single(x)*beta_q4; res_q4_alt=single(y)-yhat_q4_alt; disp('beta_q4-beta8') beta_q4-beta8 newx8_q=x; newx4_q=single(x); newx8_q(:,7)=res_q8; newx4_q(:,7)=res_q4; c8_q=corr(newx8_q); c4_q=corr(newx4_q); disp('alt c4_q via beta') newx4_q_alt=single(x); newx4_q_alt(:,7)=res_q4_alt; c4_q_alt=corr(newx4_q_alt); testres=[res8 res_q8 res4 res_q4 res_q4_alt]; testaccr =[c8(7,:)' c8_q(7,:)' c4(7,:)' c4_q(7,:)' c4_q_alt(7,:)']'; testaccr disp('++++++++++++++++ Matlab Ending +++++++++++++++') % Next two commands pass Matlab data back to Matrix % makeb34s('e3.dat',e3); % makeb34s('evec3.dat',evec3); quit % testres=[res8 res_q8 res4 res_q4] % Next two commands pass Matlab data back to Matrix % makeb34s('e3.dat',e3); % makeb34s('evec3.dat',evec3); % quit b34sreturn; %b34sendif; b34srun; == ==TLOGIT Tests Logit Model b34sexec options ginclude('b34sdata.mac') macro(murder)$ b34seend$ B34SEXEC PROBIT $ MODEL D1 = T Y LF NW $ B34SRUN$ B34SEXEC MARS LOGIT MI=2 nmcv=2 ntcv=2 $ MODEL D1=T Y LF NW $ b34srun$ b34sexec matrix; call loaddata; call load(dispmars :staging); call load(tlogit :staging); call echooff; call mars(d1 t y lf nw :nmcv 2 :ntcv 2 :mi 2 :logit :print); call names; call dispmars; call tabulate(%y %yhat); upper=.5; lower=.5; iprint=1; call character(cc,'Tests on MARS Model'); call tlogit(%y,%yhat,upper,lower,cc,ntruer,ntruep nfalser,nfalsep,nunclear,ptruer,pfalser,iprint); upper=grid(.5,.9,.01); n=norows(upper); ptruer =array(n:); pfalser=array(n:); do i=1,n; call tlogit(%y,%yhat,upper(i),lower,cc,ntruer,ntruep nfalser,nfalsep,nunclear,pt,pf,iprint); ptruer(i) =pt; pfalser(i)=pf; enddo; call graph(ptruer upper :plottype xyplot); b34srun; == ==TRIM Remove obs from fronf of a series b34sexec matrix; call load(trim :staging); call echooff; call print('Using a non time series':); x=dfloat(integers(24)); n=norows(x); newx=x; call trim(n-3,newx); call tabulate(x,newx); call graph(newx); call print('Using a time series':); x=array(:1 2 3 4 5 6 7 8 9 10); x=dfloat(integers(24)); call settime(x,1947,1,12.); call describe(x :print); dis=chardatemy(makejul(x)); call tabulate(dis,x); n=norows(x); newx=x; call trim(n-13,newx); call describe(newx :print); dis=chardatemy(makejul(x)); call tabulate(dis,x); dis=chardatemy(makejul(newx)); call tabulate(dis,newx); call print(x,newx); call graph(newx :plottype timeplot); b34srun; == ==TS_TESTS Further Tests on Residuals of OLS Model %b34slet dosur =0; %b34slet dorats=0; b34sexec options ginclude('b34sdata.mac') member(grunfeld); b34srun; b34sexec matrix; call loaddata; call load(b_g_test); call load(lmtest); call load(ts_tests :staging); call load(reset69); call echooff; /; /; Set parameters for tests /; nacf=4; bg_order=4; ireset1=1; ireset2=2; ireset3=3; lmorder=4; call print(' ':); call print('GE Equation':); call print('-----------':); call print(' ':); call olsq(i_ge f_ge c_ge :diag :print :savex); call print(' ':); call describe(%res :print); call ts_tests(nacf,%res,%y,%x,bg_order,ireset1,ireset2,ireset3,lmorder); /; call alt_aic(%nob,%k,%rss); /; do i=1,bg_order; /; call b_g_alt(i,%x,%res); /; enddo; call print(' ':); call print('Westinghouse Equation':); call print('---------------------':); call print(' ':); call olsq(i_west f_west c_west :diag :print :savex); call print(' ':); call describe(%res :print); call ts_tests(nacf,%res,%y,%x,bg_order,ireset1,ireset2,ireset3,lmorder); /; call alt_aic(%nob,%k,%rss); /; do i=1,bg_order; /; call b_g_alt(i,%x,%res); /; enddo; b34srun; %b34sif(&dosur.ne.0)%then; b34sexec simeq printsys reduced ols ls3 kcov=diag $ heading=('Theil (1970) pages 294-302' ) $ exogenous constant f_ge c_ge f_west c_west ; endogenous i_ge i_west; model lvar=i_ge rvar=(constant f_ge c_ge) name=('GE') $ model lvar=i_west rvar=(constant f_west c_west) name=('West.') $ b34seend $ %b34sendif; %b34sif(&dorats.ne.0)%then; b34sexec options open('rats.dat') unit(28) disp=unknown$ b34srun$ b34sexec options open('rats.in') unit(29) disp=unknown$ b34srun$ b34sexec options clean(28)$ b34srun$ b34sexec options clean(29)$ b34srun$ b34sexec pgmcall$ rats passasts pcomments('* ', '* Data passed from B34S(r) system to RATS', '* ', "display @1 %dateandtime() @33 ' Rats Version ' %ratsversion()" '* ') $ PGMCARDS$ * * Rats Manual User's Guide Version 5 page 190 Breusch-Godfrey (1978) Test * Rats approach is not consistent with Greene which suggests that lag * residual values that are NA be set = 0. Rats drops these obs * linreg i_ge / resids # constant f_ge c_ge stats resids linreg resids # constant f_ge c_ge resids{1} cdf chisqr %trsq 1 b34sreturn$ b34srun $ b34sexec options close(28)$ b34srun$ b34sexec options close(29)$ b34srun$ b34sexec options /$ dodos(' rats386 rats.in rats.out ') dodos('start /w /r rats32s rats.in /run') dounix('rats rats.in rats.out')$ B34SRUN$ b34sexec options npageout WRITEOUT('Output from RATS',' ',' ') COPYFOUT('rats.out') dodos('ERASE rats.in','ERASE rats.out','ERASE rats.dat') dounix('rm rats.in','rm rats.out','rm rats.dat') $ B34SRUN$ %b34sendif; == ==TSALIGN Align Time Series Data /; /; Shows line up and purging time series data. /; Due to possible missing data inside the series the timestart /; and timebase have not been set. However a date variable has been /; added to preserve the date of each observation as has been done /; here /; b34sexec matrix; call echooff; call tsd(:load :file 'c:\b34slm\examples\tsd3.tsd' /; :print :nodate :nomessage); call load(tsalign :staging); call tsalign; call tabulate(cdate,year,c,cd); call graph(year c :plottype xyplot); b34srun; == ==TSMOOTH Illustrates smooth command /; /; Illustrates "Automatic Methods" /; Results can be graphed /; /; Uses /; /; subroutine tsmooth(x,minlag,maxlag,ismooth, /; iprintdt,igraph,iprint,sum_tab,igrid); /; /; using defaults look at performance of "Automatic Methods" /; in in-sample "forecasting" /; /; x - series /; minlag - minumum lag considered /; maxlag - maximum lag considered /; ismooth - = 0 use Hankie models /; = 1 use Gardner models /; iprintdt - ne 0 print %actual and %xhat / %s /; igraph - ne 0 => draw graphs. Use minlag=maxlag /; iprint - ne 0 show smooth output /; sum_tab - ne 0 print a summary table /; igrid - ne 0 show a grid in graphs. /; /; Built 21 July 2011 /; b34sexec options ginclude('gas.b34'); b34srun; b34sexec matrix; call loaddata; call load(tsmooth :staging); call echooff; x=gasout; /; Tests for lag=1-4 forecasts minlag=1; maxlag=4; iprintdt=0; igraph=0; iprint=0; igrid=1; sum_tab=1; ismooth=0; call tsmooth(x,minlag,maxlag,ismooth,iprintdt,igraph,iprint, sum_tab,igrid); ismooth=1; call tsmooth(x,minlag,maxlag,ismooth,iprintdt,igraph,iprint, sum_tab,igrid); b34srun; b34sexec options ginclude('b34sdata.mac') member(bj_g); b34srun; b34sexec matrix; call loaddata; call load(tsmooth :staging); call echooff; x=airline; /; Tests for lag=1 forecasts minlag=1; maxlag=4; iprintdt=0; igraph=0; iprint=0; igrid=1; sum_tab=1; ismooth=0; call tsmooth(x,minlag,maxlag,ismooth,iprintdt,igraph,iprint, sum_tab,igrid); ismooth=1; call tsmooth(x,minlag,maxlag,ismooth,iprintdt,igraph,iprint, sum_tab,igrid); b34srun; == ==UNDIF Undifference a series b34sexec matrix; call load(undif :staging); call echooff; c=array(:10 20 40 80 160 320 640 1280); d_1_1=dif(c); d_1_2=dif(c,1,2); d_2_1=dif(c,2,1); /; first dif i=integers(norows(c)-1); test_1_1=c(i+1)-c(i); /; second order dif i=integers(norows(c)-2); test_1_2=c(i+2)-c(i); /; two first dif i=integers(norows(c)-2); test_2_1=test_1_1(i+1)-test_1_1(i); call tabulate(c,d_1_1,d_1_2,d_2_1,test_1_1,test_1_2,test_2_1); /; subroutine undif(dseries,order,base,nseries); /; /; Will undifference /; /; dseries => differenced series /; order => order /; base => base value of original series. Must be length of order /; nseries => new series /; call undif(d_1_1,1,c(1),new_1_1); /; now look at 2_1 base=d_1_1(1); call undif(d_2_1,1,base ,hold); base=c(1); call undif(hold, 1,base ,new2_1); /; /; looking at 1_2 /; base=c(integers(2)); call undif(d_1_2,2,base,new1_2); call tabulate(c,new_1_1,new2_1,new1_2); b34srun; == ==UPMATRIX Unpack a matrix /; /; Application of LTS to MARS and GAM Models /; b34sexec options ginclude('gas.b34'); b34srun; b34sexec matrix; call loaddata; call load(lts :staging); call load(upmatrix :staging); call load(catname :staging); call echooff; n=6; p=.9; iprint=1; call olsq(gasout gasin{1 to n} gasout{1 to n} :print :savex); call lts(%y,%x,%res,%names,%lag,%coef,%se,%t,p,newy,newx,ihold, iprint); call marspline(newy newx :print :mi 2 :nk 20 :holdout ihold); call upmatrix(newx,%names,%lag,level(),newnames); add ='[predictor,3]'; call catname(newnames,add,augname,1); call gamfit(newy argument(augname) :holdout ihold :print); b34srun; == ==WALD Tests Restrictions /; /; Problem From Greene (4th Ed) page 276 /; b34sexec matrix; call load(wald :staging); y=vector(:.161 .172 .158 .173 .195 .217 .199 .163 .195 .231 .257 .259 .225 .241 .204); x=matrix(15,5: 1. 1. 1.058 5.16 4.40 1. 2. 1.088 5.87 5.15 1. 3. 1.086 5.95 5.37 1. 4. 1.122 4.88 4.99 1. 5. 1.186 4.50 4.16 1. 6. 1.254 6.44 5.75 1. 7. 1.246 7.83 8.82 1. 8. 1.232 6.25 9.31 1. 9. 1.298 5.50 5.21 1. 10. 1.370 5.46 5.83 1. 11. 1.439 7.46 7.40 1. 12. 1.479 10.28 8.64 1. 13. 1.474 11.77 9.31 1. 14. 1.503 13.42 9.44 1. 15. 1.475 11.02 5.99); /; call olsq(y x :noint :qr :print); /; r=matrix(3,5:0 1 0 0 0 0 0 1 0 0 0 0 0 1 1); q=vector(:0 1 0); /; /; want to test b(2) = 0 /; b(3) = 1 /; b(4)+b(5)=0 /; iprint=1; call wald(%coef,%xpxinv,r,q,%nob,%resvar,waldtest, probwald,iprint); b34srun; ==