! ****************************************************************** ! staging2.mac ! Copyright(c) 2000 ! ! Created: 7/17/2009 9:24:19 PM ! Author : HOUSTON STOKES ! Last change: HS 12/17/2012 4:00:33 AM ! ****************************************************************** ==README /; /; Library of subroutines and functions that extend the capability /; of the Matrix Command. These commands are not documented in the /; b34s help file for the matrix command. The commands in staging2.mac /; have internal documentation. /; /; The command /; /; call load(ace_ols :staging); /; will load the ace_ols program that is called as: /; /; call ace_ols; /; /; The ace_ols program calles the internal subroutine aceols2 /; that could be called directly. /; == ==BOOST BOOST, BOOST2, BOOST3 & BOOST4 code /; /; These routines are subject to change /; /; BOOST does the usual Boosting procedure /; BOOST2 does Modified Forward Statewise model boosting /; BOOST3 Works simular to BOOST but does not require data centering /; saves vector info and coefficients for using with boost4 /; BOOST4 Allows out of sample forecasting for model estimated with /; boost3 /; subroutine boost(y,yhat,res,x,in,e,ipass,itype,iprint); /; /; 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 /; /; y => left hand variable /; yhat => Forecasted y /; res => error /; x => n by k matrix of right hand side variables /; in => work vector telling what vectors in X used /; in = 0 implies that vector in /; e => Step size - use a small number /; ipass => Number of times called. Be sure called enough times /; one way to test this is to monitor progress in the /; correlation improvement /; itype => =0 OLS, =1 MARS, =2 GAM, =3 L1 =4 mimimax /; iprint => =0 no printing, =1 print step data, =2 print correlations /; /; Note: X must be centered. Y must have mean removed. Code to do this /; is shown next /; /; 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; /; e=.01 ; /; ntry=900; /; /; itype=0 => ols /; itype=1 => mars /; itype=2 => gam not ready /; itype=0; /; iprint=0; /; /; x=center2(catcol(age sex bmi bp s1 s2 s3 s4 s5 s6)); /; y=y-mean(y); /; /; 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); /; 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; /; /; /; +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /; This routine only adds one tree at a time. /; /; This procedure is very effective is data mining applications /; where there are many columns in X and the objective is to /; obtain the information in x without the "cost" of forming /; X'X /; /; +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /; cc=array(nocols(x):); if(ipass.eq.1)then; in=array(nocols(x):)+1.; yhat=y*0.0; res=y; endif; do i=1,nocols(x); cc(i)=ccf(res,x(,i)); enddo; /; cc=afam(cc)*afam(in); ij=imax(abs(cc)); if(iprint.eq.2)then; call print('Largest Correlation vector was ',ij:); call tabulate(in,cc); endif; /; in(ij)=0.0; if(iprint.eq.0)then; if(itype.eq.0)call olsq( res,x(,ij) :noint); if(itype.eq.1)call marspline(res,x(,ij) ); if(itype.eq.2)call gamfit( res,x(,ij)[predictor,3] :noint); if(itype.eq.3)call olsq( res,x(,ij) :noint :l1); if(itype.eq.4)call olsq( res,x(,ij) :noint :minimax); endif; if(iprint.ne.0)then; if(itype.eq.0)call olsq( res,x(,ij) :noint :print); if(itype.eq.1)call marspline(res,x(,ij) :print ); if(itype.eq.2)call gamfit( res,x(,ij)[predictor,3] :noint :print); if(itype.eq.3)call olsq( res,x(,ij) :noint :print :l1); if(itype.eq.4)call olsq( res,x(,ij) :noint :print :minimax); endif; if(itype.eq.3)%yhat=%l1yhat; if(itype.eq.4)%yhat=%mmyhat; yhat=yhat+ afam(e)*afam(%yhat); res=y-yhat; return; end; subroutine boost2(y,yhat,res,x,xbuild,in,e,ipass,itype,iprint); /; /; Implements modified OLS/MARS/GAM boosting /; Suggested by below listed reference as an interesting extension /; to boosting. This can be thought of as Modified Forward Stagewise /; boosting /; /; Modified OLS boosting is described 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 /; /; y => left hand variable /; yhat => Forecasted y /; res => error /; x => n by k matrix of right hand side variables /; in => work vector telling what vectors in X used in set = 0 /; e => Step size /; ipass => Number of times called /; itype => =0 OLS, =1 MARS, =2 GAM, =3 L1, =4 Minimax /; iprint => =0 no printing, =1 print step data, =2 print correlations /; /; Note: X must be centered. Y must have mean removed. Code to do this /; is shown next /; /; 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; /; e=.01 ; /; ntry=900; /; /; itype=0 => ols /; itype=1 => mars /; itype=2 => gam not ready /; itype=0; /; iprint=0; /; /; x=center2(catcol(age sex bmi bp s1 s2 s3 s4 s5 s6)); /; y=y-mean(y); /; /; 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); /; 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; /; /; /; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /; This routine only adds new tree plus all old trees at each step. /; This is called modified Forward Statewise /; /; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /; cc=array(nocols(x):); if(ipass.eq.1)then; in=array(nocols(x):)+1.; yhat=y*0.0; res=y; endif; do i=1,nocols(x); cc(i)=ccf(res,x(,i)); enddo; /; cc=afam(cc)*afam(in); ij=imax(abs(cc)); if(iprint.eq.2)then; call print('Largest Correlation vector was ',ij:); call tabulate(in,cc); endif; if(iprint.ne.0)call print(ipass,cc,ij); if(ipass.eq.1) xbuild=array(norows(x),1:x(,ij)); if(ipass.gt.1.and.in(ij).ne.0.0)xbuild=catcol(xbuild,x(,ij)); in(ij)=0.0; do j=1,nocols(xbuild); if(iprint.eq.0)then; if(itype.eq.0)call olsq( res,xbuild(,j):noint); if(itype.eq.1)call marspline(res,xbuild(,j)); if(itype.eq.2)call gamfit( res,xbuild(,j)[predictor,3] :noint); if(itype.eq.3)call olsq( res,xbuild(,j):noint :l1); if(itype.eq.4)call olsq( res,xbuild(,j):noint :minimax); endif; if(iprint.ne.0)then; if(itype.eq.0)call olsq( res,xbuild(,j) :noint :print); if(itype.eq.1)call marspline(res,xbuild(,j) :print); if(itype.eq.2)call gamfit( res,xbuild(,j)[predictor,3] :noint :print); if(itype.eq.3)call olsq( res,xbuild(,j) :noint :print :l1); if(itype.eq.4)call olsq( res,xbuild(,j) :noint :print :minimax); endif; if(itype.eq.3)%yhat=%l1yhat; if(itype.eq.4)%yhat=%mmyhat; yhat=yhat+ afam(e)*afam(%yhat); res=y-yhat; enddo; return; end; subroutine boost3(y,yhat,res,x,in,beta1,beta2,e,ipass,itype,iprint); /; /; Implements a mod to OLS boosting as described 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 /; /; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /; HHS Mods involve use of constant and no mean adjustment /; Goal is to facilitate forecasting using boost4 routine /; in, beta1, beta2 used as input into boost4 /; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /; /; y => left hand variable /; yhat => Forecasted y /; res => error /; x => n by k matrix of right hand side variables /; in => work vector telling vector in X used /; beta1 => saves the beta /; beta2 => sames the constant /; e => Step size - use a small number /; ipass => Number of times called. Be sure called enough times /; one way to test this is to monitor progress in the /; correlation improvement /; itype => =0 OLS, =1 MARS, =2 GAM, =3 L1 =4 mimimax /; iprint => =0 no printing, =1 print step data, =2 print correlations /; /; b34sexec matrix; /; call loaddata; /; call load(center :staging); /; call load(center2 :staging); /; call load(boost :staging); /; call echooff; /; /; e=.01 ; /; ntry=900; /; /; itype=0 => ols /; itype=1 => mars not ready /; itype=2 => gam not ready /; itype=3 => L1 /; itype=4 => MINIMAX /; itype=0; /; iprint=0; /; /; 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 outdouble(20,3,fit(i)); /; enddo; /; /; /; +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /; This routine only adds one tree at a time. /; /; This procedure is very effective is data mining applications /; where there are many columns in X and the objective is to /; obtain the information in x without the "cost" of forming /; X'X /; /; +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /; cc=array(nocols(x):); if(itype.eq.1.or.itype.eq.2)then; call print('ERROR: ITYPE in boots3 cannot be 1 or 2':); call stop; endif; if(ipass.eq.1)then; yhat=y*0.0; res=y; endif; do i=1,nocols(x); cc(i)=ccf(res,x(,i)); enddo; ij=imax(abs(cc)); in(ipass)=ij; if(iprint.eq.2)then; call print('Largest Correlation vector was ',ij:); call tabulate(in,cc); endif; if(iprint.eq.0)then; if(itype.eq.0)call olsq( res,x(,ij) ); /; if(itype.eq.1)call marspline(res,x(,ij) ); /; if(itype.eq.2)call gamfit( res,x(,ij)[predictor,3] :noint); if(itype.eq.3)call olsq( res,x(,ij) :l1); if(itype.eq.4)call olsq( res,x(,ij) :minimax); endif; if(iprint.ne.0)then; if(itype.eq.0)call olsq( res,x(,ij) :print); /; if(itype.eq.1)call marspline(res,x(,ij) :print ); /; if(itype.eq.2)call gamfit( res,x(,ij)[predictor,3] :noint :print); if(itype.eq.3)call olsq( res,x(,ij) :print :l1); if(itype.eq.4)call olsq( res,x(,ij) :print :minimax); endif; if(itype.eq.0)then; beta1(ipass)=%coef(1); beta2(ipass)=%coef(2); endif; if(itype.eq.3)then; %yhat=%l1yhat; beta1(ipass)=%l1coef(1); beta2(ipass)=%l1coef(2); endif; if(itype.eq.4)then; %yhat=%mmyhat; beta1(ipass)=%mmcoef(1); beta2(ipass)=%mmcoef(2); endif; yhat=yhat+ afam(e)*afam(%yhat); res=y-yhat; return; end; subroutine boost4(yhat,x,in,beta1,beta2,e); /; /; HHS boosting Forecasting for model estimated with boost3 /; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /; HHS Mods involve use of constant and no mean adjustment /; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /; /; yhat => Forecasted y /; x => n by k matrix of right hand side variables /; in => work vector telling what vectors in X used /; beta1 => Beta Coef /; beta2 => Constant coef /; e => Step size used /; /; +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /; yhat=array(norows(x):); do i=1,norows(beta1); add=(beta1(i)*x(,in(i)))+beta2(i); yhat=yhat+ afam(e)*afam(add); enddo; return; end; == ==CATNAME Add string to a names vector subroutine catname(left,right,newname,drop); /; /; Augments the name vector to add [predictor,3] etc /; /; Usually run on a namelist which is character*8 and has 1 or more /; elements. /; /; Can be used to add a suffex to a character*1 array /; /; left => character*8 name list or character*1 scaler /; right => string to put at end /; newname => New augmented name vector/array /; drop => # elements to drop from front for 1d left values /; /; Built 21 June 2007 by Houston H. Stokes /; Mod 7 January 2008 that reversed arguments and allowed subsets of /; oldname /; /; ------------------------------------------------------------- /; if(kind(left).eq.-1.or.(kind(left).eq.-8.and.klass(left).eq.0))then; call character(left2,left); call character(right2,right); newname=catrow(left2,right2); go to done; endif; if(kind(left).eq.-8.and.klass(left).ne.0)then; if(klass(left).ne.5.and.klass(left).ne.1)then; call epprint('ERROR: CATNAME does not support 2D input arrays':); go to done; endif; if(drop.ge.norows(left))go to done; call character(right2,right); i1 =norows(right2); i2 =8+1; newname=c1array(norows(left)-drop,i1+i2:); jj =integers(1,i1+i2); do i=1,norows(left)-drop; call character(left2,left(i)); fill=catrow(left2,right2); newname(i,)=fill; enddo; go to done; endif; done continue; if(klass(newname).eq.6.and.nocols(newname).eq.1) newname=to_array(newname); if(klass(newname).eq.2.and.nocols(newname).eq.1) newname=to_vector(newname); if(klass(newname).eq.2.or.klass(newname).eq.6) newname=transpose(newname); return; end; == ==CENTER Center a series or 2 D object mean=0 variance=1 function center(x); /; /; Standardizes data such that mean = 0 and unit variance /; X can be an array, vector or matrix /; /; Built 24 February 2008 by Houston H. Stokes if(kind(x).ne.8.and.kind(x).ne.-16)then; call epprint('ERROR: function center requires real*8 or real*16 data':); dataout=missing(); go to finish; endif; if(norows(x).le.1)then; call epprint('ERROR: function center passed one element object':); dataout=missing(); go to finish; endif; dataout=afam(x); if(nocols(x).eq.1)then; v=variance(dataout); dataout=dataout-mean(dataout); if(v.ne.kindas(x,0.0))dataout=dataout/sqrt(v); if(klass(x).eq.1)dataout=vfam(dataout); endif; if(nocols(x).gt.1)then; do i=1,nocols(x); v=variance(dataout(,i)); dataout(,i)=(dataout(,i)-mean(dataout(,i))); if(v.ne.kindas(x,0.0)) dataout(,i)=(dataout(,i))/sqrt(v); enddo; if(klass(x).eq.2)dataout=mfam(dataout); endif; finish continue; return(dataout); end; == ==CENTER2 Standardizes a Series mean=0 Unit length function center2(x); /; /; Standardizes data such that mean = 0 and series has unit length /; X can be an array, vector or matrix /; /; Built 24 February 2008 by Houston H. Stokes if(kind(x).ne.8.and.kind(x).ne.-16)then; call epprint('ERROR: function center2 requires real*8 or real*16 data': ); dataout=missing(); go to finish; endif; if(norows(x).le.1)then; call epprint('ERROR: function center2 passed one element object':); dataout=missing(); go to finish; endif; dataout=afam(x); if(nocols(x).eq.1)then; dataout=dataout-mean(dataout); v=sumsq(dataout); if(v.ne.kindas(x,0.0))dataout=(dataout-mean(dataout))/sqrt(v); if(klass(x).eq.1)dataout=vfam(dataout); endif; if(nocols(x).gt.1)then; do i=1,nocols(x); dataout(,i)=(dataout(,i)-mean(dataout(,i))); v=sumsq(dataout(,i)); if(v.ne.kindas(x,0.0))dataout(,i)=(dataout(,i))/sqrt(v); enddo; if(klass(x).eq.2)dataout=mfam(dataout); endif; finish continue; return(dataout); end; == ==CGETR16 Reads Real*16 Data into Matrix subroutine cgetr16(n_row,n_col,x16); /; /; Reads data directly into real*16 array(n_row,n_col) /; data on unit 4 /; /; n_row = # of rows /; n_col = # of cols /; x16 = data read /; call rewind(4); x16=r8tor16(array(n_col,n_row:)); call read(x16,4); x16=transpose(x16); call rewind(4); return; end; == ==CGETVPA Reads VPA data into Matrix subroutine cgetvpa(n_row,n_col,xvpa); /; /; Reads data directly into vpa array(n_row,n_col) /; data on unit 4 /; /; n_row = # of rows /; n_col = # of cols /; xvpa = data read /; /; Requires call load(getvpa :staging); /; call rewind(4); cc=c1array(10000:); call read(cc,4); nn=n_row*n_col; call getvpa(cc,nn,xvpa,ibad); xvpa=array(n_row,n_col:xvpa); call rewind(4); return; end; == ==COR Correlation of a matrix. function cor(x); /; /; Use matrix language to correlate y1 and y2 which are cols x /; Can use real*4, real*8, real*16 and VPA. For real*8 use ccf(x,y) /; /; Need to call load(cov :staging); /; k_hold=klass(x); if(norows(x).le.nocols(x).or.(klass(x).ne.6.and.klass(x).ne.2))then; call epprint('ERROR: Matrix passed to cor not correct':); call epprint(' # rows ',norows(x):); call epprint(' # cols ',nocols(x):); call epprint(' kind ',klass(x):); call stop; endif; ccov=afam(cov(x)); n1=norows(ccov); n2=nocols(ccov); ccor=array(n1,n2:); zero=kindas(x,0.0); do i=1,n1; do j=1,n2; if(ccov(i,i).ne.zero.and.ccov(j,j).ne.zero) ccor(i,j)=ccov(i,j)/sqrt(ccov(i,i)*ccov(j,j)); enddo; enddo; return(ccor); end; == ==COV Calculated Covariance of a matrix function cov(x); /; /; Use matrix language to calculate cov of a matrix. For series use /; mm=catcol(series1,series2) /; Can use real*4, real*8, real*16 and VPA. /; test=afam(x); i=nocols(x); d=kindas(x,dfloat(norows(x)-1)); do j=1,i; test(,j)=test(,j)-mean(test(,j)); enddo; ccov=afam(transpose(mfam(test))*mfam(test))/d; if(klass(x).eq.2.or.klass(x).eq.1)ccov=mfam(ccov); return(ccov); end; == ==COV2 Idempotent Matrix method get Variance-Cov function cov2(x); /; /; Use matrix language to calculate cov of a matrix. For series use /; mm=catcol(series1,series2) /; Can use real*4, real*8, real*16 and VPA. /; /; Uses Greene(2000,16) idempotent M0 matrix /; function cov( ) is a more traditional approach that uses /; far less space at the cost of a speed loss due to do loops /; cov( ) is more accurate than cov2( ) if there are scaling /; differences /; /; At issue is that m0 is n by n !!!!! /; /; Use of i which is a vector of 1's /; z=catcol(x1,x2,...,xk) /; ccov=transpose(z)*m0*z/(1/(n-1)); /; where m0 diagonal = 1-(1/n). Off Diag = -1/n /; n=norows(x); real_n=kindas(x,dfloat(n)); real_one=kindas(x,1.0); /; Define i as a n by 1 matrix of ones i=matrix(n,1:kindas(x,(vector(n:)+1.))); /; Get Variance Covariance /; Define Idempotent matrix M Diagonal = 1-(1/n). Off Diag -(1/n) bigi=kindas(x,matrix(n,n:)) + real_one; littlei=(real_one/real_n)*(i*transpose(i)); m0=(bigi-littlei); ccov2=transpose(mfam(x))*m0*mfam(x)/(real_n-real_one); if(klass(x).eq.6.or.klass(x).eq.5)ccov2=afam(ccov2); return(ccov2); end; == ==COND Condition of a Matrix function cond(x,i); /; Condition of a matrix /; /; This command tracks Matlab 100% /; /; cond(x,i) = norm(x,i)*norm(inv(x),i) /; /; x = square matrix; /; i = code in range 1-4 /; cond(x,1) largest norm1 col sum /; cond(x,2) largest singular value of x / smallest singular value /; cond(x,3) largest norm1 row sum /; cond(x,4) Frobenius norm sqrt(sum(diag(x'*x))) /; /; Built 4 April 2007 /; /; tracks matlab 100% /; /; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ if(kind(i).ne.-4)then; call epprint('ERROR: cond(x,i) must pass i as an integer'); call stop; endif; if(klass(x).ne.2)then; call epprint('ERROR: cond( ) must pass matrix as argument # 1'); call stop; endif; if(i.lt.1.or.i.gt.4)then; call epprint('ERROR: cond(a,i) for matrix a must set i=1,...,4'); call stop; endif; if(i.ne.2)then; ans=norm(x,i)*norm(inv(x),i); endif; if(i.eq.2)then; k=nocols(x); s=svd(x); if(s(k).eq.kindas(x,0.0))then; ans=kindas(x,.1e+32); endif; if(s(k).ne.kindas(x,0.0))then; ans=s(1)/s(k); endif; endif; return(ans); end; == ==DDOT dot product of vectors with an offset 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 /; /; Note: dot(norows(x),x,1,x,1) & prod(x) give same answer /; 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; == ==DF_PP Dickey Fuller & Phillips Perron Tests subroutine df_pp(x,info,n); /; /; Battery of DF, PP, ERS and KPSS tests /; /; x => series to test /; title => Title /; n => # of lags /; Routine built 19 October 2008 /; ERS added April 2009 /; /; Note that ppunit and dfunit 100% track Hamilton and SAS /; KPSS tracks Greene and Rats /; need /; call load(kpss); /; call load(df_gls); /; /; This command gives t and z forms of DF and PP. The command /; /; call uroot_t(x,title,n,iprint,makeg); /; /; gives a battery of DF, PP, ERS and KPSS tests for a unit root /; /; x => series to test /; title => Title /; n => # of lags /; iprint => 0 no printing. Used with makeg=1 /; 1 print table /; 2 detail as well as table /; -1 same as 1 except critical value table not printed /; -2 same as 2 except critical value table not printed /; makeg => 1 makes global %lag %df_i, %adf, %adft,%pp_i,%app, /; %appt %kpssnt,%kpsst,%ers_nt, %ers_t, /; %df1%, %df5%, %df10%, /; %adf1%, %adf5%, %adf10%, /; %adft1%,%adft5%,%adft10%); /; /; /; uroot_t does not give the z form of the statistics but it is more /; compact and gives critical values. /; /; --------------------------------------------------------- pp_I =array(n+1:); df_I =array(n+1:); pp_Iz =array(n+1:); df_Iz =array(n+1:); app =array(n+1:); appt =array(n+1:); lag=array(n+1:); app_z =array(n+1:); appt_z =array(n+1:); adf =array(n+1:); adft =array(n+1:); adf_z =array(n+1:); adft_z =array(n+1:); kpssnt=array(n+1:); kpsst =array(n+1:); ers_nt=array(n+1:); ers_t =array(n+1:); do i=0,n; j=i+1; call df(x, a1 :df i); df_i(j)=a1; call df(x, a1 :df i :zform); df_iz(j)=a1; call df(x, a1 :adf i); adf(j) =a1; call df(x, a1 :adf i :zform); adf_z(j) =a1; call df(x, a2 :adft i); adft(j) =a2; call df(x, a2 :adft i :zform); adft_z(j)=a2; call pp(x, a1 :pp i); pp_i(j)=a1; call pp(x, a1 :pp i :zform); pp_iz(j)=a1; call pp(x, a1 :app i); app(j) =a1; call pp(x, a1 :app i :zform); app_z(j) =a1; call pp(x, a2 :appt i); appt(j) =a2; call pp(x, a2 :appt i :zform); appt_z(j)=a2; call kpss(x,test1,test2,i,0); kpssnt(j)=test1; kpsst(j) =test2; call df_gls(x,i,_notrend,_trend,notrendx,trendx,0); ers_t(j) =_trend; ers_nt(j) =_notrend; lag(j)=dfloat(i); enddo; call print('DF, Phillips-Perron, KPSS & ERS tests on ',title:); mm=mean(x); vv=variance(x); ss=sqrt(vv); nn=norows(x); call print('Mean of Series ',mm:); call print('Variance of Series ',vv:); call print('Stardard Deviation of Series ',ss:); call print('Number of Observations in Series',nn:); call print('df_i => df ':); call print('df_iz => df z form ':); call print('adf => augmented df ':); call print('adf_z => augmented df ':); call print('adft => augmented df with trend':); call print('adft_z=> augmented df with trend':); call print('pp_i => pp test ':); call print('pp_iz => pp test z form ':); call print('app => augmented pp ':); call print('app_z => augmented pp ':); call print('appt => augmented pp with trend':); call print('appt_z=> augmented pp with trend':); call print(' ':); call print( 'KPSS Stationarity Test without trend (kpssnt) Critical Velues':); call print('Critical Values .10 .050 .025 .010 ':); call print(' .347 .463 .573 .739 ':); call print(' ':); call print( 'KPSS Stationarity Test with trend (kpsst) Critical Values':); call print('Critical Values .10 .050 .025 .010 ':); call print(' .119 .146 .176 .216 ':); call print(' ':); call print('ers_nt critical values are 10% -1.62, 5% -1.95, 1% -2.58':); call print('ers_t critical values are 10% -2.57, 5% -2.89, 1% -3.48':); call print('Test value < critical value => reject unit root ':); call print(' ':); call tabulate(lag,df_i,df_iz,adf,adf_z,adft,adft_z,kpssnt,kpsst); call tabulate(lag,pp_i,pp_iz,app,app_z,appt,appt_z,ers_nt,ers_t); return; end; == ==DFVALUES Estimate of DF/PP Critical Values subroutine dfvalues(iopt,cvalue,noob); /; /; Obtains an estimate of the DF Critical Values /; designed to track RATS 100% /; iopt =1 => Intercept and trend /; iopt =2 => Intercept /; iopt =3 => Neither Intercept or trend /; values(1) => 1% /; values(2) => 2% /; values(3) => 3% /; noob => # of obs /; critical values of Dickey Fuller /; Fuller, Introduction to Statistical Time Series, New York, Wiley, /; 1976. Dickey and Fuller, "Distribution of the Estimators for Time /; Series Regressions with a Unit Root", J.A.S.A., 1979, pp 427-431. /; The (approximate) critical values for t-test form are from /; MacKinnon, "Critical Values for Cointegration Tests", Long-Run /; Economic Relationships, R.F. Engle and C.W.J. Granger, eds, /; London, Oxford, 1991, pp 267-276 /; /; Based on Internal B34S Subroutine based on Rats Logic /; /; Built 10 March 2007 /; --------------------------------------------------------------- xnoob=dfloat(noob); cvalue=array(3:); /; /; With intercept and trend /; if(iopt.eq.1)then; cvalue(1)=-3.9638 - (8.353/xnoob) - 47.44/(xnoob**2.); cvalue(2)=-3.4126 - (4.039/xnoob) - 17.83/(xnoob**2.); cvalue(3)=-3.1279 - (2.418/xnoob) - 7.58/(xnoob**2.); endif; /; /; With intercept /; if(iopt.eq.2)then; cvalue(1)=-3.4335 - (5.999/xnoob)-(29.25/(xnoob**2.)); cvalue(2)=-2.8621 - (2.738/xnoob)-( 8.36/(xnoob**2.)); cvalue(3)=-2.5671 - (1.438/xnoob)-( 4.48/(xnoob**2.)); endif; /; /; With neither intercept or trend /; if(iopt.eq.3)then; cvalue(1)=-2.5658 -(1.960/xnoob)- (10.04/(xnoob**2.)); cvalue(2)=-1.9393 -(0.398/xnoob); cvalue(3)=-1.6156 -(0.181/xnoob); endif; /; return; end; == ==DISPMARS Displays Mars Model program dispmars; /; /; Displays MARS Output - Used after call mars( ); /; Modified Nov 5, 2005 (BL) /; Modified December 2005 by Houston H. Stokes to run only if /; MARS (1991) has been called /; Changed 16 January by Houston H. Stokes to save i /; /; Displays both the MARS model estimated and the relative importance /; of variables in MARS model. /; /; Relative Importance Section /; save_i=0; if(kind(i).ne.-99)then; hold_i=i; save_i=1; endif; save_j=0; if(kind(j).ne.-99)then; hold_j=j; save_j=1; endif; save_kk=0; if(kind(kk).ne.-99)then; hold_kk=kk; save_kk=1; endif; save_kkk=0; if(kind(kkk).ne.-99)then; hold_kkk=kkk; save_kkk=1; endif; save_nn=0; if(kind(nn).ne.-99)then; hold_nn=nn; save_nn=1; endif; if(norows(%names).gt.1.and.%mars_vr.eq.0)then; gvsort=array(norows(%names):); gvsort =ranker(%varrimp); gvindx = integers(1,norows(%names)); gvindx = norows(%names)-gvindx+1; varimp=%varrimp(gvsort); vnames=%names(gvsort); vtypes=%typevar(gvsort); minval=%minvar(gvsort); maxval=%maxvar(gvsort); vlag=%lag(gvsort); varimp=varimp(gvindx); vnames=vnames(gvindx); vtypes=vtypes(gvindx); minval=minval(gvindx); maxval=maxval(gvindx); vlag=vlag(gvindx); fmtmin=c1array(norows(%names):); do i=1,norows(%names); if(dabs(minval(i)).ge.1.0)then; if(dabs(minval(i)).lt.1000.00000000)fmtmin(i)='(f18.8)'; if(dabs(minval(i)).ge.1000.00000000)fmtmin(i)='(f18.8)'; if(dabs(minval(i)).ge.10000.0000000)fmtmin(i)='(f18.7)'; if(dabs(minval(i)).ge.100000.000000)fmtmin(i)='(f18.6)'; if(dabs(minval(i)).ge.1000000.00000)fmtmin(i)='(f18.5)'; if(dabs(minval(i)).ge.10000000.0000)fmtmin(i)='(f18.4)'; if(dabs(minval(i)).ge.100000000.000)fmtmin(i)='(f18.3)'; if(dabs(minval(i)).ge.1000000000.00)fmtmin(i)='(f18.2)'; if(dabs(minval(i)).ge.10000000000.0)fmtmin(i)='(f18.1)'; if(dabs(minval(i)).gt.10000000000.0)fmtmin(i)='(g18.2)'; endif; if(dabs(minval(i)).lt.1.0)then; if(dabs(minval(i)).lt.0.000000000001)fmtmin(i)='(g18.6)'; if(dabs(minval(i)).ge.0.000000000001)fmtmin(i)='(f18.12)'; if(dabs(minval(i)).ge.0.000000000010)fmtmin(i)='(f18.11)'; if(dabs(minval(i)).ge.0.000000000100)fmtmin(i)='(f18.10)'; if(dabs(minval(i)).ge.0.000000001000)fmtmin(i)='(f18.9)'; if(dabs(minval(i)).ge.0.000000010000)fmtmin(i)='(f18.8)'; if(dabs(minval(i)).eq.0.000000000000)fmtmin(i)='(f18.8)'; endif; enddo; fmtmax=c1array(norows(%names):); do i=1,norows(%names); if(dabs(maxval(i)).ge.1.0)then; if(dabs(maxval(i)).lt.1000.00000000)fmtmax(i)='(f18.8)'; if(dabs(maxval(i)).ge.1000.00000000)fmtmax(i)='(f18.8)'; if(dabs(maxval(i)).ge.10000.0000000)fmtmax(i)='(f18.7)'; if(dabs(maxval(i)).ge.100000.000000)fmtmax(i)='(f18.6)'; if(dabs(maxval(i)).ge.1000000.00000)fmtmax(i)='(f18.5)'; if(dabs(maxval(i)).ge.10000000.0000)fmtmax(i)='(f18.4)'; if(dabs(maxval(i)).ge.100000000.000)fmtmax(i)='(f18.3)'; if(dabs(maxval(i)).ge.1000000000.00)fmtmax(i)='(f18.2)'; if(dabs(maxval(i)).ge.10000000000.0)fmtmax(i)='(f18.1)'; if(dabs(maxval(i)).gt.10000000000.0)fmtmax(i)='(g18.2)'; endif; if(dabs(maxval(i)).lt.1.0)then; if(dabs(maxval(i)).lt.0.000000000001)fmtmax(i)='(g18.6)'; if(dabs(maxval(i)).ge.0.000000000001)fmtmax(i)='(f18.12)'; if(dabs(maxval(i)).ge.0.000000000010)fmtmax(i)='(f18.11)'; if(dabs(maxval(i)).ge.0.000000000100)fmtmax(i)='(f18.10)'; if(dabs(maxval(i)).ge.0.000000001000)fmtmax(i)='(f18.9)'; if(dabs(maxval(i)).ge.0.000000010000)fmtmax(i)='(f18.8)'; if(dabs(maxval(i)).eq.0.000000000000)fmtmax(i)='(f18.8)'; endif; enddo; call fprint(:clear :col 1 :string ' ' :print); call fprint(:clear :col 1 :string 'Relative Importance of Variables in MARS model' :print); call fprint(:clear :col 1 :string '----------------------------------------------' :print); call fprint(:clear :col 1 :display 'Variable' :col 14 :display 'Importance' :col 25 :display 'VarType' :col 38 :display 'Lag' :col 54 :display 'Minimum' :col 74 :display 'Maximum' :print ); do i=1,norows(vnames); if(varimp(i).le.0.)exitdo; if(vtypes(i).eq.0)call character(vtmp,'Continuous '); if(vtypes(i).ne.0)call character(vtmp,'Categorical'); call fprint(:clear :col 1 :display vnames(i) :col 12 :display varimp(i) '(f8.2)' :col 26 :display vtmp :col 38 :display vlag(i) :col 44 :display minval(i) fmtmin :col 64 :display maxval(i) fmtmax :print); enddo; endif; /; Works only for MARS(1991) if(%mars_vr.eq.0)then; call fprint(:clear :col 1 :cr 1 :string 'Mars Results after Stepwise Elimination' :print :cr 1 :clear :col 1 :string %yvar :col 10 :string '=' :col 16 :display %coef(1) '(g16.8)' :print :clear); do j=2,norows(%coef); if(%coef(j).ne.0.0)then; /; ********************************************************* /; coef*max(var - value,0.0)) /; coef*max(value - var,0.0) /; if(%typevar(%varink(j-1)).eq.0)then; call fprint(:clear :col 1 :display %coef(j) '(g16.8)' ); if(%coef(j).ge.0.0)call fprint(:col 1 :string '+'); ibase=18; call disp_max(j-1,%typek,%names,%varink,%lag,%knot, %typevar,%states,ibase,0.0); call fprint(:print); /; if(%parent(j-1).ne.0)then; /; kk=%parent(j-1); /; ibase=18; /; call fprint(:clear); /; call disp_max(kk,%typek,%names,%varink,%lag,%knot, /; %typevar,%states,ibase,0.0); /; call fprint(:print); /; if(%parent(kk).ne.0)then; /; kkk=%parent(kk); /; ibase=18; /; call fprint(:clear); /; call disp_max(kkk,%typek,%names,%varink,%lag,%knot, /; %typevar,%states,ibase,0.0); /; call fprint(:print); /; endif; /; endif; if(%parent(j-1).ne.0)then; kk=%parent(j-1); if(%typevar(%varink(kk)).ne.0)then; chold2=%cknot(kk,); call ialen(chold2,jjchold); nn=0; ibase=18; do nntest=jjchold,1,-1; nn=nn+1; call character(one,'1'); if(chold2(nntest).eq.one(1))then; call fprint(:clear); ibase=18; call disp_max(kk,%typek,%names,%varink,%lag,%knot, %typevar,%states,ibase,dfloat(nntest)); call fprint(:print); endif; enddo; endif; if(%typevar(%varink(kk)).eq.0)then; call fprint(:clear); ibase=18; call disp_max(kk,%typek,%names,%varink,%lag,%knot, %typevar,%states,ibase,0.0); call fprint(:print); endif; if(%parent(kk).ne.0)then; kkk=%parent(kk); if(%typevar(%varink(kkk)).ne.0)then; chold3=%cknot(kkk,); call ialen(chold3,jjchold); nn=0; ibase=18; do nntest=jjchold,1,-1; nn=nn+1; call character(one,'1'); if(chold3(nntest).eq.one(1))then; call fprint(:clear); ibase=18; call disp_max(kkk,%typek,%names,%varink,%lag,%knot, %typevar,%states,ibase,dfloat(nntest)); call fprint(:print); endif; enddo; endif; if(%typevar(%varink(kkk)).eq.0)then; call fprint(:clear); ibase=18; call disp_max(kkk,%typek,%names,%varink,%lag,%knot, %typevar,%states,ibase,0.0); call fprint(:print); endif; endif; endif; endif; /; ******************************************************** /; if(var .eq. value) /; if(%typevar(%varink(j-1)).ne.0)then; chold=%cknot(j-1,); call ialen(chold,jjchold); nn=0; do nntest=jjchold,1,-1; nn=nn+1; call character(one,'1'); if(chold(nntest).eq.one(1))then; call fprint(:clear :col 1 :display %coef(j) '(g16.8)' ); if(%coef(j).ge.0.0)call fprint(:col 1 :string '+'); ibase=18; call disp_max(j-1,%typek,%names,%varink,%lag,%knot, %typevar,%states,ibase,dfloat(nntest)); call fprint(:print); endif; enddo; if(%parent(j-1).ne.0)then; kk=%parent(j-1); if(%typevar(%varink(kk)).ne.0)then; chold2=%cknot(kk,); call ialen(chold2,jjchold); nn=0; ibase=18; do nntest=jjchold,1,-1; nn=nn+1; call character(one,'1'); if(chold2(nntest).eq.one(1))then; call fprint(:clear); ibase=18; call disp_max(kk,%typek,%names,%varink,%lag,%knot, %typevar,%states,ibase,dfloat(nntest)); call fprint(:print); endif; enddo; endif; if(%typevar(%varink(kk)).eq.0)then; call fprint(:clear); ibase=18; call disp_max(kk,%typek,%names,%varink,%lag,%knot, %typevar,%states,ibase,0.0); call fprint(:print); endif; if(%parent(kk).ne.0)then; kkk=%parent(kk); if(%typevar(%varink(kkk)).ne.0)then; chold3=%cknot(kkk,); call ialen(chold3,jjchold); nn=0; ibase=18; do nntest=jjchold,1,-1; nn=nn+1; call character(one,'1'); if(chold3(nntest).eq.one(1))then; call fprint(:clear); ibase=18; call disp_max(kkk,%typek,%names,%varink,%lag,%knot, %typevar,%states,ibase,dfloat(nntest)); call fprint(:print); endif; enddo; endif; if(%typevar(%varink(kkk)).eq.0)then; call fprint(:clear); ibase=18; call disp_max(kkk,%typek,%names,%varink,%lag,%knot, %typevar,%states,ibase,0.0); call fprint(:print); endif; endif; endif; endif; /; ******************************************************** endif; enddo; endif; if(save_i.eq.1) i=hold_i; if(save_j.eq.1) j=hold_j; if(save_kk.eq.1) kk=hold_kk; if(save_kkk.eq.1)kkk=hold_kkk; if(save_nn.eq.1) nn=hold_nn; return; end; subroutine disp_max(jj,itype,names,varink,lag,knot,typevar, states,ibase,cknoti); /; /; jj => # for all but coef which is plus 1 /; itype => %typek 0 => var-knot, 1 => knot-var /; names => %names names vector /; varink => Variable number /; lag => Lag of variable /; knot => knot value for continuous case. State # for 0-1 /; case except if more than 2 states /; typevar see below /; typevar=0 => Displays coef itype 0 * max(var-knot,0.0) /; typevar=0 => Displays coef itype 1 * max(knot-var,0.0) /; typevar NE 0 => Displays coef itype 0 if(var eq knot) /; typevar ne 0 => Displays coef itype 1 if(var ne knot) /; states => State values for 0-1 variables. missing otherwise /; ibase => Offset to write /; cknoti => states # from character*1 array of 28 elements /; showing values that are on for 0-1 case /; /; Program revised 20 October 2004 /; ********************************************************************* /; max(var-knot,0.0) vame=c1array(8:names(varink(jj))); call ialen(vame,ii); if(ii.lt.8)then; vame=moveright(vame,(8-ii)); endif; if(itype(jj).eq.0.and.typevar(varink(jj)).eq.0)then; call fprint(:col ibase :string '* max(' :col 24 :string vame :col 32 :string '{' :col 33 :display lag(varink(jj)) '(i3)' :col 36 :string '} -' :col 40 :display knot(jj) '(g16.8)' :col 56 :string ' , 0.0)'); endif; if(itype(jj).eq.0.and.typevar(varink(jj)).ne.0)then; value=missing(); call getcat(typevar,states,cknoti,varink(jj),value); if(value.ne.missing())then; call fprint(:col ibase :string 'if(' :col 24 :string vame :col 32 :string '{' :col 33 :display lag(varink(jj)) '(i3)' :col 36 :string '} eq' :col 40 :display value '(g16.8)' :col 56 :string ')'); endif; if(value.eq.missing())then; call fprint(:col ibase :string 'if (' :col 24 :string vame :col 32 :string '{' :col 33 :display lag(varink(jj)) '(i3)' :col 36 :string '} eq' :col 40 :string ' na value' :col 56 :string ')'); endif; endif; /; max(knot-var,0.0) if(itype(jj).eq.1.and.typevar(varink(jj)).eq.0)then; call fprint(:col ibase :string '* max(' :col ibase+6 :display knot(jj) '(g16.8)' :col ibase+22 :string ' - ' :col ibase+26 :string vame :col ibase+34 :string '{' :col ibase+36 :display lag(varink(jj)) '(i3)' :col ibase+39 :string '}, 0.0)'); endif; if(itype(jj).eq.1.and.typevar(varink(jj)).ne.0)then; value=missing(); call getcat(typevar,states,cknoti,varink(jj),value); if(value.ne.missing())then; call fprint(:col ibase :string 'if(' :col ibase+6 :display value '(g16.8)' :col ibase+22 :string ' eq' :col ibase+26 :string vame :col ibase+34 :string '{' :col ibase+36 :display lag(varink(jj)) '(i3)' :col ibase+39 :string '})'); endif; if(value.eq.missing())then; call fprint(:col ibase :string 'if(' :col ibase+6 :string ' NA values' :col ibase+22 :string ' - ' :string ' eq' :col ibase+26 :string vame :col ibase+34 :string '{' :col ibase+36 :display lag(varink(jj)) '(i3)' :col ibase+39 :string '}'); endif; endif; return; end; subroutine getcat(typevar,states,knot,varinki,value); /; /; Get Cat variable value /; /; typevar see below /; typevar=0 => Displays coef itype 0 * max(var-knot,0.0) /; typevar=0 => Displays coef itype 1 * max(knot-var,0.0) /; typevar NE 0 => Displays coef itype 0 if(var eq knot) /; typevar ne 0 => Displays coef itype 1 if(var ne knot) /; states => State values for 0-1 variables. missing otherwise /; knot => knot value for continuous case. State # for 0-1 /; case except if more than 2 states /; varinki => variable in knot /; /; Program revised 20 October 2004 /; ***************************************************************** isum=0; if(varinki.gt.1)then; do i=1,varinki-1; isum=isum+typevar(i); enddo; endif; value=missing(); ii=isum+idint(knot); if(ii.gt.0)then; value=states(isum+idint(knot)); endif; return; end; == ==F_MARSLG Adjusts yhat and forecasts from MARSPLINE subroutine f_marslg(raw_s,adj1,adj2); /; /; Normalizes a prediction in the range 0-1 /; /; raw_s => the output from matspline /; adj1 => Values > 1.0 and < 0.0 set to 1.0 and 0.0 /; adj2 => Series rescaled /; /; This routine is 100% experimental /; +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ adj1 =dmin1(raw_s,1.0); adj1 =dmax1(adj1 ,0.0); mm1=min(raw_s); mm2=max(raw_s); range=mm2-mm1; if(mm1.lt.0.0)then; range=mm2+dabs(mm1); adj2 =raw_s+dabs(mm1); endif if(mm1.gt.0.0)adj2 =raw_s -mm1; adj2 =adj2 /range; return; end; == ==G4TEST Pena-Slate JASA March (2006) Joint OLS Test /; /; g4test => g4 test driver /; g4test2( ) => calculations /; g4move => Driver for leverage plots /; g4move2( ) => calculations for leverage plots /; g4plot => plots for leverage plots /; g4rr => driver for g4 rr calculations /; g4rr2( ) => calculations for g4 rr /; g4rrplot => plots for g4 rr /; /; Full test: /; 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(%rg4dif, %rpg4dif, %rs1dif ,%rps1dif, %rs2dif, /; %rps2dif,%rs3dif, %rps3dif,%rs4dif, %rps4dif); /; b34srun; /; /; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /; program g4test; call g4test2(%res,%y,%yhat,%coef %s1,%s2,%s3,%s4,%g4,%x,%ps1,%ps2,%ps3,%ps4,%pg4,1); return; end; subroutine g4test2(res,y,yhat,beta,s1, s2, s3, s4,g4,x, ps1,ps2,ps3,ps4,pg4,iprint); /; Joint tests of OLS assumptions /; /; res => OLS Residual /; y => OLS y /; yhat => OLS yhat /; S1 => Linearity Test /; S2 => Homoscedasticity Test /; S3 => Uncorrelatedness /; S4 => Normality Test /; G4 => Joint Test /; PS1 => Probability S1 /; PS2 => Probability S2 /; PS3 => Probability S3 /; PS4 => Probability S4 /; PS5 => Probability S5 /; iprint => ser = 1 to print results /; /; Edsel Pena and Elizabeth Slate /; "Global Validation of Linear Model Assumptions" /; JASA March 2006, Vol. 101 No. 473 /; /; Built 1 June 2006 by Houston H. Stokes. /; /; Edsel Pena was most helpful in releasing R code for this problem. /; /; Examples of Use: /; /; 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; /; /; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /; n=norows(res); k=norows(beta); div=1./dfloat(n); /; small sample /; div=1./dfloat(n-nocols(x)); r = afam(res)/dsqrt(sumsq(res)*div); y = afam(y); yhat = afam(yhat); /; Stat S1 & S2 s1 =(sum( r**3.0) /dsqrt( 6.0*dfloat(norows(res))))**2.; ps1 =chisqprob(s1,1.); s2 =(sum((r**4.0)-3.0)/dsqrt(24.0*dfloat(norows(res))))**2.; ps2 =chisqprob(s2,1.); /; Stat S3 yh_meany=yhat-mean(y); s3 =((1.0/dsqrt(dfloat(n)))*sum(yh_meany*yh_meany*r))**2.; b1=(sum((yh_meany)**4.))*div; newx=afam(array(n,k-1:)); do i=1,k-1; newx(,i)=x(,i)-mean(x(,i)); enddo; b22= transpose(mfam(newx))*mfam(newx)*div; newbeta=beta; newbeta(k)=missing(); newbeta=goodrow(newbeta); b2=vfam(newbeta)*mfam(b22)*vfam(newbeta); b3=mfam((afam(yh_meany)*afam(yh_meany)*afam(newx))*div); gama=vector(k-1:); do i=1,k-1; gama(i)=sum(b3(,i)); enddo; b3=gama*inv(b22)*gama; s3=s3/(b1-(b2*b2)-b3); ps3 =chisqprob(s3,1.); v=dfloat(integers(n))/dfloat(n); mean_v=mean(v); /; stat S4 s4 =(sum((v-mean_v)*((r*r)-1.0))/ dsqrt(2.0*sum((v-mean_v)*(v-mean_v))))**2.; ps4 =chisqprob(s4,1.); g4 =s1+s2+s3+s4; pg4 =chisqprob(g4,4.); if(iprint.ne.0)then; call print(' ':); call print('Pena - Slate (JASA March 2006) Specification Tests:':); call print( 'Test Statistic Probability':); call fprint(:clear :col 1 :string 'S1 - Linearity ' :col 30 :display s1 '(g12.4)' :col 50 :display ps1 '(g12.4)' :print); call fprint(:clear :col 1 :string 'S2 - Homoscedasticity ' :col 30 :display s2 '(g12.4)' :col 50 :display ps2 '(g12.4)' :print); call fprint(:clear :col 1 :string 'S3 - Uncorrelatedness ' :col 30 :display s3 '(g12.4)' :col 50 :display ps3 '(g12.4)' :print); call fprint(:clear :col 1 :string 'S4 - Normality Test ' :col 30 :display s4 '(g12.4)' :col 50 :display ps4 '(g12.4)' :print); call fprint(:clear :col 1 :string 'G4 - Joint Test ' :col 30 :display g4 '(g12.4)' :col 50 :display pg4 '(g12.4)' :print); endif; return; end; subroutine g4move_2(g4base,pg4base, s1base,ps1base, s2base,ps2base, s3base,ps3base, s4base,ps4base, g4dif, pg4dif, s1dif, ps1dif, s2dif, ps2dif, s3dif, ps3dif, s4dif, ps4dif, y,x); /; /; Calculates Outlier Detection Statistics suggested by Pena-Slate /; /; Edsel Pena and Elizabeth Slate /; "Global Validation of Linear Model Assumptions" /; JASA March 2006, Vol. 101 No. 473 /; /; Built 1 June 2006 by Houston H. Stokes. /; /; Edsel Pena was most helpful in releasing R code for this problem. /; /; /; g4base => g4 statistic for complete sample /; pg4base => probability g4 statistic for complete sample /; g4dif =array(norows(x):)+g4base; pg4dif=array(norows(x):)+pg4base; s1dif =array(norows(x):)+s1base; ps1dif=array(norows(x):)+ps1base; s2dif =array(norows(x):)+s2base; ps2dif=array(norows(x):)+ps2base; s3dif =array(norows(x):)+s3base; ps3dif=array(norows(x):)+ps3base; s4dif =array(norows(x):)+s4base; ps4dif=array(norows(x):)+ps4base; do i=1,norows(x); xnew=x; ynew=y; xnew(i,)=missing(); ynew(i) =missing(); xnew=goodrow(xnew); ynew=goodrow(ynew); call olsq(ynew xnew :noint :savex); iprint=0; call g4test2(%res,%y,%yhat,%coef,s1, s2, s3, s4,g4,%x, ps1,ps2,ps3,ps4,pg4,iprint); g4dif(i) =(g4-g4dif(i))/g4dif(i); s1dif(i) =(s1-s1dif(i))/s1dif(i); s2dif(i) =(s2-s2dif(i))/s2dif(i); s3dif(i) =(s3-s3dif(i))/s3dif(i); s4dif(i) =(s4-s4dif(i))/s4dif(i); pg4dif(i)= pg4; ps1dif(i)= ps1; ps2dif(i)= ps2; ps3dif(i)= ps3; ps4dif(i)= ps4; enddo; g4dif=g4dif*100.; s1dif=s1dif*100.; s2dif=s2dif*100.; s3dif=s3dif*100.; s4dif=s4dif*100.; return; end; subroutine g4rr_2(rg4, rpg4,rs1,rps1,rs2,rps2,rs3,rps3,rs4,rps4,y,x); /; /; Calculates Outlier Detection Statistics suggested by Pena-Slate /; Uses recursive calculation /; /; Edsel Pena and Elizabeth Slate /; "Global Validation of Linear Model Assumptions" /; JASA March 2006, Vol. 101 No. 473 /; /; Built 1 June 2006 by Houston H. Stokes. /; /; Edsel Pena was most helpful in releasing R code for this problem. /; /; n=norows(x); k=nocols(x); rg4 =array(n-k-1:); rpg4=array(n-k-1:); rs1 =array(n-k-1:); rps1=array(n-k-1:); rs2 =array(n-k-1:); rps2=array(n-k-1:); rs3 =array(n-k-1:); rps3=array(n-k-1:); rs4 =array(n-k-1:); rps4=array(n-k-1:); icount=0; do i=k+2,n; xnew=x; ynew=y; icount=icount+1; j=integers(i,n); xnew(j,1)=missing(); ynew(j) =missing(); xnew=goodrow(xnew); ynew=goodrow(ynew); call olsq(ynew xnew :noint :savex); iprint=0; call g4test2(%res,%y,%yhat,%coef,s1, s2, s3, s4, g4,%x, ps1, ps2, ps3, ps4,pg4,iprint); rg4(icount) =g4; rs1(icount) =s1; rs2(icount) =s2; rs3(icount) =s3; rs4(icount) =s4; rpg4(icount)=pg4; rps1(icount)=ps1; rps2(icount)=ps2; rps3(icount)=ps3; rps4(icount)=ps4; enddo; return; end; program g4move; /; /; Drive g4move_2 /; /; Must have :savex on /; call olsq( ) /; /; and have a prior call /; call g4test; /; /; Example: /; 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); y=%y; x=%x; call g4move_2(%g4,%pg4,%s1,%ps1,%s2,%ps2,%s3,%ps3,%s4,%ps4, %g4dif, %pg4dif, %s1dif ,%ps1dif, %s2dif, %ps2dif, %s3dif, %ps3dif, %s4dif, %ps4dif, y,x); return; end; program g4rr; /; /; Drive g4rr_2 /; /; Must have :savex on /; call olsq( ) /; /; and have a prior call /; call g4test; /; /; Example: /; call olsq(n_gal miles_lf daysbetw :print :savex ); /; call g4test; /; call g4rr; /; call g4rrplot; /; call tabulate(%rg4, %rpg4, %rs1 , %rps1, %rs2, /; %rps2,%rs3, %rps3, %rs4, %rps4); y=%y; x=%x; call g4rr_2(%rg4, %rpg4, %rs1 ,%rps1, %rs2, %rps2, %rs3, %rps3, %rs4, %rps4, y,x); return; end; program g4plot; /; /; Plot g4move data /; call graph(%g4dif :heading 'Change in G4 Test ' :xlabel 'Obs Dropped' 'C' :pgborder :pgborder :pgyscaleleft 'NT' :pgyscaleright 'I' :pgxscaletop 'I' :nocontact :file 'g4dif.wmf'); call graph(%pg4dif :heading 'Change in G4 Test Probability ' :xlabel 'Obs Dropped' 'C' :pgborder :pgyscaleleft 'NT' :pgyscaleright 'I' :pgxscaletop 'I' :nocontact :file 'pg4dif.wmf'); call graph(%s1dif :heading 'Change in S1 Test (Linearity) ' :xlabel 'Obs Dropped' 'C' :pgborder :pgyscaleleft 'NT' :pgyscaleright 'I' :pgxscaletop 'I' :nocontact :file 's1dif.wmf'); call graph(%ps1dif :heading 'Change in S1 Test Probability (Linearity)' :xlabel 'Obs Dropped' 'C' :pgborder :pgyscaleleft 'NT' :pgyscaleright 'I' :pgxscaletop 'I' :nocontact :file 'ps1dif.wmf'); call graph(%s2dif :heading 'Change in S2 Test (Homoscedasticity)' :xlabel 'Obs Dropped' 'C' :pgborder :pgyscaleleft 'NT' :pgyscaleright 'I' :pgxscaletop 'I' :nocontact :file 's2dif.wmf'); call graph(%ps2dif :heading 'Change in S2 Test Probability (Homoscedasticity)' :xlabel 'Obs Dropped' 'C' :pgborder :pgyscaleleft 'NT' :pgyscaleright 'I' :pgxscaletop 'I' :nocontact :file 'ps2dif.wmf'); call graph(%s3dif :heading 'Change in S3 Test (Uncorrelatedness)' :xlabel 'Obs Dropped' 'C' :pgborder :pgyscaleleft 'NT' :pgyscaleright 'I' :pgxscaletop 'I' :nocontact :file 's3dif.wmf'); call graph(%ps3dif :heading 'Change in S3 Test Probability (Uncorrelatedness) ' :xlabel 'Obs Dropped' 'C' :pgborder :pgyscaleleft 'NT' :pgyscaleright 'I' :pgxscaletop 'I' :nocontact :file 'ps3dif.wmf'); call graph(%s4dif :heading 'Change in S4 Test (Normality)' :xlabel 'Obs Dropped' 'C' :pgborder :pgyscaleleft 'NT' :pgyscaleright 'I' :pgxscaletop 'I' :nocontact :file 's4dif.wmf'); call graph(%ps4dif :heading 'Change in S4 Test Probability (Normality) ' :xlabel 'Obs Dropped' 'C' :pgborder :pgyscaleleft 'NT' :pgyscaleright 'I' :pgxscaletop 'I' :nocontact :file 'ps4dif.wmf'); return; end; program g4rrplot; /; /; Plot g4rr data /; call graph(%rg4 :heading 'Recursive G4 Test ' :xlabel 'Obs Added Recursively' 'C' :pgborder :pgyscaleleft 'NT' :pgyscaleright 'I' :pgxscaletop 'I' :nocontact :file 'rg4dif.wmf'); call graph(%rpg4 :heading 'Recursive G4 Test Probability ' :xlabel 'Obs Added Recursively' 'C' :pgborder :pgyscaleleft 'NT' :pgyscaleright 'I' :pgxscaletop 'I' :nocontact :file 'rpg4dif.wmf'); call graph(%rs1 :heading 'Recursive S1 Test (Linearity)' :xlabel 'Obs Added Recursively' 'C' :pgborder :pgyscaleleft 'NT' :pgyscaleright 'I' :pgxscaletop 'I' :nocontact :file 'rs1dif.wmf'); call graph(%rps1 :heading 'Recursive S1 Test Probability (Linearity) ' :xlabel 'Obs Added Recursively' 'C' :pgborder :pgyscaleleft 'NT' :pgyscaleright 'I' :pgxscaletop 'I' :nocontact :file 'rps1dif.wmf'); call graph(%rs2 :heading 'Recursive S2 Test (Homoscedasticity)' :xlabel 'Obs Added Recursively' 'C' :pgborder :pgyscaleleft 'NT' :pgyscaleright 'I' :pgxscaletop 'I' :nocontact :file 'rs2dif.wmf'); call graph(%rps2 :heading 'Recursive S2 Test Probability (Homoscedasticity)' :xlabel 'Obs Added Recursively' 'C' :pgborder :pgyscaleleft 'NT' :pgyscaleright 'I' :pgxscaletop 'I' :nocontact :file 'rps2dif.wmf'); call graph(%rs3 :heading 'Recursive S3 Test (Uncorrelatedness) ' :xlabel 'Obs Added Recursively' 'C' :pgborder :pgyscaleleft 'NT' :pgyscaleright 'I' :pgxscaletop 'I' :nocontact :file 'rs3dif.wmf'); call graph(%rps3 :heading 'Recursive S3 Test Probability (Uncorrelatedness)' :xlabel 'Obs Added Recursively' 'C' :pgborder :pgyscaleleft 'NT' :pgyscaleright 'I' :pgxscaletop 'I' :nocontact :file 'rps3dif.wmf'); call graph(%rs4 :heading 'Recursive S4 Test (Normality)' :xlabel 'Obs Added Recursively' 'C' :pgborder :pgyscaleleft 'NT' :pgyscaleright 'I' :pgxscaletop 'I' :nocontact :file 'rs4dif.wmf'); call graph(%rps4 :heading 'Recursive S4 Test Probability (Normality) ' :xlabel 'Obs Added Recursively' 'C' :pgborder :pgyscaleleft 'NT' :pgyscaleright 'I' :pgxscaletop 'I' :nocontact :file 'rps4dif.wmf'); return; end; == ==GARCHF Forecast ARCH / GARCH Series Subroutine garchf(mtype,res1,res2,y,nfore,fore1,fore2,iobs, nparm,cname,coef,corder,nar,nma,ngar,ngma,nmu,con, nosqrt,iprintf); /; /; Note: This routine is experimental and subject to change /; -------------------------------------------------- /; /; Generate Forecasts /; /; Arguments: /; mtype => Model type not used in release # 1 /; res1 => first moment residual /; res2 => second moment residual /; y => series /; nfore => # forecasts /; fore1 => first moment forecast /; fore2 => second moment forecast /; iobs => Observations of the forecast /; nparm => Number of parameters in model /; cname => Internal names of the model paraemters /; coef => ar, ma, gar,gma,my, vd, const1,const1 /; corder => Coef orders /; nar => # ar terms /; nma => # ma terms /; ngar => # second moment ar terms /; ngma => # second monent ma terma /; nmu => # mu terms /; nosqrt => =0 Sqrt in GARCH-M; =1 NoSqrt in GARCH-M /; con => # constants - coded -1, 0, 1 2 /; iprintf=> =1 print, =2 fprint, =0 no print /; /; Note: Usually nar nma ngar ngma con set by garchest as /; %nar %nma %ngar %nnma %con /; /; Preliminary version May 30, 2005 by Houston H. Stokes /; use with caution /; /; Modified on July 26, 2005 (BL) /; + IGARCH & GARCH-M logic added /; + Exit if regression variables are in model (not supported) /; + Added logic to allow for non-normal distributions for forecasting /; /; Forecast logic from "Analysis of Financial Time Series" by Ruey Tsay /; /; /; Note: ARMA, ARCH, GARCH, IGARCH, and GARCH-M supported at this time. /; /; arma uses Tsay page 53 f_y(1) = const1+ sum ar + sum MA /; arch uses Tsay page 90 f_res2(1)=a0+a1*res2(n)+a2*res2(n-1) /; f_res2(2)=a0+a1*f_res2(1)+a2*res2(n) /; garch(1,1) uses Tsay page 94 f_res2(1)=a0 /; /; fore1=array(nfore+1:); fore2=array(nfore+1:); iobs= array(nfore+1:integers(norows(res1),norows(res1)+nfore)); /; /; Logic: we call garch but pass extended data from forecast. /; /; call garch(res2,arch2,g,func,2,n /; :ar array(:b1,b2) idint(array(:1 2)) /; :ma array(: /; :gma array(:a1) idint(array(:1) ) /; :gar array(:a2) idint(array(:1) ) /; :constant array(:b0 a0) /; :forecast norows(g) ); arg1=' '; arg2=' '; arg3=' '; arg4=' '; arg5=' '; arg6=' '; arg7=' '; icount=0; /; If there are input variables - EXIT (NOT SUPPORTED)!!! strinp='INPUT'; strtmp=c1array(8:strinp); do i=1,nparm; strvnm=c1array(8:cname(i)); ifind=0; do j=1,5; if(strtmp(j).eq.strvnm(j)) ifind=ifind+1; enddo; if(ifind.eq.5) icount=icount+1; if(ifind.eq.5)then; call print('** Regression+GARCH Forecasting not supported':); return; endif; enddo; if(nar .ne.0)i1=integers(icount+1,icount+nar); icount=icount+nar; if(nma .ne.0)i2=integers(icount+1,icount+nma); icount=icount+nma; if(ngar.ne.0)i3=integers(icount+1,icount+ngar); icount=icount+ngar; if(ngma.ne.0)i4=integers(icount+1,icount+ngma); icount=icount+ngma; if(nmu.ne.0)i5=integers(icount+1,icount+nmu); icount=icount+nmu; /; Count the number of VD terms (fattail) /; Forecasts are the same as Normal distribution do i=1,nparm; if(cname(i).eq.'VD') icount=icount+1; enddo; if(con .eq.2) i6=integers(icount+1,icount+con); if( abs(con) .eq.1)i6=integers(icount+1,icount+1); icount=icount+abs(con); if(nar .ne.0) arg1=':ar array(:coef(i1)) idint(array(:corder(i1)))'; if(nma .ne.0) arg2=':ma array(:coef(i2)) idint(array(:corder(i2)))'; if(ngar.ne.0) arg3=':gar array(:coef(i3)) idint(array(:corder(i3)))'; if(ngma.ne.0) arg4=':gma array(:coef(i4)) idint(array(:corder(i4)))'; if(nmu .ne.0) arg5=':mu array(:coef(i5)) idint(array(:corder(i5)))'; if(con .eq.2) arg6=':constant array(:coef(i6)) '; if(con .eq.-1)arg6=':constant array(:coef(i6),0.0) '; if(con .eq.1) arg6=':constant array(:0.0,coef(i6))'; if(nosqrt.eq.1)arg7=':nosqrt'; maxlag=1; if(nar .ne.0)maxlag=max1(max(idint(array(:corder(i1)))),maxlag); if(nma .ne.0)maxlag=max1(max(idint(array(:corder(i2)))),maxlag); if(ngar.ne.0)maxlag=max1(max(idint(array(:corder(i3)))),maxlag); if(ngma.ne.0)maxlag=max1(max(idint(array(:corder(i4)))),maxlag); if(nmu .ne.0)maxlag=max1(max(idint(array(:corder(i5)))),maxlag); if(con .ne.0)maxlag=max1(max(idint(array(:corder(i6)))),maxlag); gg=y; rres1=res1; rres2=res2; if(mtype.eq.3)then; /; IGARCH model is being forecasted - build constrained gma1=1-gar1 pgmav=1.-coef(i3); arg4=':gma array(:pgmav) idint(array(:corder(i3)))'; endif; do ii=1,nfore+1; /; call testarg(rres1,rres2,gg,func,maxlag,n /; Note: Must not have ; after argument( ) since these are options!! /; call garch(rres1,rres2,gg,func,maxlag,n argument(arg1) argument(arg2) argument(arg3) argument(arg4) argument(arg5) argument(arg6) argument(arg7) :forecast norows(y)); n=norows(gg)+1; gg(n) =%f_m1_m2(1); rres2(n)=%f_m1_m2(2); fore1(ii)=%f_m1_m2(1); fore2(ii)=%f_m1_m2(2); enddo; if(iprintf.eq.1)call tabulate(iobs,fore1,fore2); if(iprintf.eq.2)then; call prntfcst(norows(res1),fore1,fore2); endif; return; end; subroutine prntfcst(iobs,fore1,fore2); /$ /$ Prints GARCH Forecast Table using FPRINT routine /$ IOBS (INPUT) => Number of observations for residuals /$ FORE1 (INPUT) => Array of First moment forecasts /$ FORE2 (INPUT) => Array of Second moment forecasts /$ nrows=norows(fore1); iobsy=iobs; call print(' ':); * call fprint(:col 1 * :string 'Forecast Table for First Moment and Second Moment' * :print :clear); call print( 'Forecast Table for First Moment and Second Moment':); call print( '-----------------------------------------------------':); /$ 123456789012345678901234567890123456789012345678901234567890 call print( 'FCST Y-INDEX LEVEL FCST VOLATILITY FCST':); call print( '----- ---------- ---------------- ----------------':); *call fprint(:col 1 :string '----------' * :col 10 :string '----------' * :col 20 :string '----------' * :col 30 :string '----------' * :col 40 :string '--------------' * :print :clear); *call fprint(:col 1 :string 'NFCST' * :col 8 :string ' Y-INDEX' * :col 20 :string ' LEVEL FCST' * :col 38 :string ' VOLATILITY FCST' * :print :clear); *call fprint(:col 1 :string '-----' * :col 8 :string '----------' * :col 20 :string '----------------' * :col 38 :string '----------------' * :print :clear); jj=0; do i=1,nrows; call fprint(:col 1 :display jj '(I5)' :col 8 :display iobsy '(I10)' :col 20 :display fore1(i) '(f16.4)' :col 38 :display fore2(i) '(f16.5)' :print :clear); jj=jj+1; iobsy=iobsy+1; enddo; call print(' ':); return; end; == ==GARCH2PA Automatic GARCH Two Pass Modeling 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 /; if(noprint.ne.0)then; nt1=nt; if(nt1.le.0) nt1=norows(x); if(nt1.gt.norows(x))nt1=norows(x); nf1=nf; if(nf.le.0)nf1=1; if(nf1.gt.50)nf1.50; call autobj(x,nac,res1 :autobuild :print :nac nac :npac nac :seasonal seas :forecast index(nf1,nt1) ); res1 =%res; ressq1=%res*%res; /; nt2=norows(ressq1); call autobj(ressq1 :autobuild :print :nac nac :npac nac :seasonal seas :forecast index(nf1,nt2) ); endif; /; if(noprint.eq.0)then; call autobj(x,nac,res1 :autobuild :nac nac :npac nac :seasonal seas ); res1 =%res; ressq1=%res*%res; /; call autobj(ressq1 :autobuild :nac nac :npac nac :seasonal seas ); endif; return; end; == ==GENGARCH Simulate GARCH and ARCH Models subroutine gengarch(n,s2,gar,gma,idrop,e,h); /; /; Routine to simulate a GARCH(p,q) or ARCH(0,q) Model /; /; n => # of observations /; s2 => Unconditional Variance /; gar => Parameters for lagged conditional Variance /; This is needed for GARCH. If ARCH is desired, /; pass 0.0 /; gma => Parameters for lagged squared innovations /; idrop => # of terms to drop /; e => Vector of GARCH(p,q) innovations /; h => Associated conditional variance process /; /; Notes: This implementation assumes that the first moment model /; does not contain ar & ma terms. /; /; Logic of routine based on GAUSS routine GARCH developed /; by Mico Loretan /; /; Version 1 August 2004 /; ******************************************************* /; if(kind(gar).ne.8.or.kind(gma).ne.8)then; call epprint('ERROR: gar and gma must be real*8'); return; endif; if(klass(gar).ne.5.or.klass(gma).ne.5)then; call epprint('ERROR: gar and gma must be 1D array types'); return; endif; if(sfam(sum(gar)+sum(gma)).gt.1.0)then; call epprint('ERROR: sum(gar)+sum(gma) > 1.0'); return; endif; if(dmin(gar).lt.0.0.or.dmin(gma).lt.0.0)then; call epprint('ERROR: element in gar or gma LT 0'); return; endif; a0=afam((1.0-sum(gar)-sum(gma))*s2); e2=rn(array(n+idrop:)) *dsqrt(s2); h2= (array(n+idrop:)+1.)*s2; gar2=vector(:vfam(gar(integers(norows(gar),1,-1)))); gma2=vector(:vfam(gma(integers(norows(gma),1,-1)))); ibegin=dmax1(norows(gar),norows(gma)); do i =ibegin+1,n+idrop; kk1=integers(i-norows(gar),i-1); kk2=integers(i-norows(gma),i-1); h22=vector(:vfam(h2(kk1))); e22=vector(:vfam(e2(kk2)*e2(kk2))); h2(i)=a0+afam(gar2*h22)+afam(gma2*e22); e2(i)=dsqrt(h2(i))*rn(1.0); enddo; j=integers(idrop+1,n+idrop); e=e2(j); h=h2(j); return; end; == ==GETCHAR Unpack a string into Character Array subroutine getchar(cc,nfind,cx,ibad); /; /; Converts a string of numbers to character*1 array /; /; cc => Character*1 string /; nfind => number of numbers in cc to load. If less found /; rest of arry set to missing /; cx => the data series saved as 2-d Character array /; ibad => 0 all ok /; 1 problem /; /; Written 1 February 2005 by Houston H. Stokes /; Useful in reading character data into an array /; /; ++++++++++++++++++++++++++++++++++++++++++++++++++ if(kind(cc).ne.-1)then; call epprint('ERROR: Argument # 1 of ntokin is not character*1') ibad=1; go to finish; endif; /; Get max size isize2=0; ibase=1; do j=1,nfind; imax=0; call ib34s11(cc,ibase,ifbase,isize,itokty,inewp,imax); ibase=inewp; /; call print('isize was ',isize); /; call print('ibase was ',ibase); /; call print('ifbase was ',ifbase); /; call print('inewp was ',inewp); if(isize.eq.0.or.inewp.eq.0)go to fin2; isize2=dmax1(isize2,isize); if(inewp.eq.-99)go to fin2; enddo; fin2 continue; cx=c1array(nfind,isize2:); ibase=1; ibad=0; do j=1,nfind; imax=0; call ib34s11(cc,ibase,ifbase,isize,itokty,inewp,imax); if(isize.eq.0)go to finish; i=integers(ifbase,ifbase+isize-1); cx(j,)=cc(i); ibase=inewp; if(inewp.eq.-99)go to finish; enddo; finish continue; return; end; == ==GETDATA Simple Real*8 data Read routine subroutine getdata(file,nobs,nvar,data); /; /; Easy to use Matrix Data Read; /; /; file => File of real*8 data with nobs rows and nvar cols /; nobs => # of cols of data /; nvar => # of nows of data /; data => Array of Data /; /; Note that the data is read by row into to a Fortran by-col /; object. /; /; assume file x1(1) x2(1) ..... xk(1) /; . /; x1(n) x2(n) ......xk(n) /; /; read as nobs=k, nvar=n then transpose /; /; Built 24 February 2007 /; /; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /; call open(71,file); call rewind(71); data=array(nobs,nvar:); call read(data,71); call rewind(71); call close(71); return; end; == ==GETDATA2 Simple Real*16 data Read routine subroutine getdata2(file,nobs,nvar,data); /; /; Easy to use Matrix Data Read; /; /; file => File of real*16 data with nobs rows and nvar cols /; nobs => # of rows of data /; nvar => # of cols of data /; data => Array of Data /; /; This reads directly into real*16 /; /; Built 24 February 2007 /; /; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /; call open(71,file); call rewind(71); data=r8tor16(array(nobs,nvar:)); call read(data,71); call rewind(71); call close(71); return; end; == ==GETR16 Unpack a string into Real*16 Array subroutine getr16(cc,nfind,x,ibad); /; /; Converts a string of numbers to real*16 array /; /; cc => Character*1 string /; nfind => number of numbers in cc to load. If less found /; rest of array set to missing /; x => the data series saved as real*16 /; ibad => 0 all ok /; 1 problem /; /; Written 1 August 2004 by Houston H. Stokes /; Useful in reading character data into an array /; /; ++++++++++++++++++++++++++++++++++++++++++++++++++ x=r8tor16(array(nfind:)+missing()); ibase=1; ibad=0; if(kind(cc).ne.-1)then; call epprint('ERROR: Argument # 1 of ntokin is not character*1') ibad=1; go to finish; endif; do j=1,nfind; imax=0; call ib34s11(cc,ibase,ifbase,isize,itokty,inewp,imax); /; call print('isize was ',isize); /; call print('ibase was ',ibase); /; call print('ifbase was ',ifbase); /; call print('inewp was ',inewp); if(isize.eq.0.or.inewp.eq.0)go to finish; if((dabs(itokty).ne.5.and.dabs(itokty).ne.6))then; call print('ERROR: Problem reading Data') call print(' Type of token found ',itokty); i=integers(ifbase,ifbase+isize-1); find=cc(i); call print(' Token found was: '); call print(find :); call print(' # of token in string was ',j); call print(' ' :); ibad=1; call epprint('ERROR: A token is not real*8 or integer') go to finish; endif; i=integers(ifbase,ifbase+isize-1); ibase=inewp; i=integers(ifbase,ifbase+isize-1); find=cc(i); x(j)=real16(find); if(inewp.eq.-99)go to finish; enddo; finish continue; return; end; == ==GETVPA Unpack a string into VPA Array subroutine getvpa(cc,nfind,x,ibad); /; /; Converts a string of numbers to VPA array /; /; cc => Character*1 string /; nfind => number of numbers in cc to load. If less found /; rest of arry set to missing /; x => the data series saved as VPA /; ibad => 0 all ok /; 1 problem /; /; Written 1 February 2005 by Houston H. Stokes /; Useful in reading character data into an array /; /; ++++++++++++++++++++++++++++++++++++++++++++++++++ x=vpa(array(nfind:)); ibase=1; ibad=0; if(kind(cc).ne.-1)then; call epprint('ERROR: Argument # 1 of ntokin is not character*1') ibad=1; go to finish; endif; do j=1,nfind; imax=0; call ib34s11(cc,ibase,ifbase,isize,itokty,inewp,imax); /; call print('isize was ',isize); /; call print('ibase was ',ibase); /; call print('ifbase was ',ifbase); /; call print('inewp was ',inewp); if(isize.eq.0.or.inewp.eq.0)go to finish; if((dabs(itokty).ne.5.and.dabs(itokty).ne.6))then; call print('ERROR: Problem reading Data') call print(' Type of token found ',itokty); i=integers(ifbase,ifbase+isize-1); find=cc(i); call print(' Token found was: '); call print(find :); call print(' # of token in string was ',j); call print(' ' :); ibad=1; call epprint('ERROR: A token is not real*8 or integer') go to finish; endif; i=integers(ifbase,ifbase+isize-1); ibase=inewp; i=integers(ifbase,ifbase+isize-1); find=cc(i); x(j)=vpa(find); if(inewp.eq.-99)go to finish; enddo; finish continue; return; end; == ==GETR8 Unpack a string into Real*8 Array subroutine getr8(cc,nfind,x,ibad); /; /; Converts a string of numbers to real*8 array /; /; cc => Character*1 string /; nfind => number of numbers in cc to load. /; If less found any excess will be set to missing /; x => the data series saved as real*8 /; ibad => 0 all ok /; 1 problem /; /; Written 1 August 2004 by Houston H. Stokes /; Useful in reading character data into an array /; /; ++++++++++++++++++++++++++++++++++++++++++++++++++ x=array(nfind:)+missing(); ibase=1; ibad=0; if(kind(cc).ne.-1)then; call epprint('ERROR: Argument # 1 of ntokin is not character*1') ibad=1; go to finish; endif; do j=1,nfind; imax=0; call ib34s11(cc,ibase,ifbase,isize,itokty,inewp,imax); /; call print('isize was ',isize); /; call print('ibase was ',ibase); /; call print('ifbase was ',ifbase); /; call print('inewp was ',inewp); if(isize.eq.0.or.inewp.eq.0)go to finish; if((dabs(itokty).ne.5.and.dabs(itokty).ne.6))then; call print('ERROR: Problem reading Data') call print(' Type of token found ',itokty); i=integers(ifbase,ifbase+isize-1); find=cc(i); call print(' Token found was: '); call print(find :); call print(' # of token in string was ',j); call print(' ' :); ibad=1; call epprint('ERROR: A token is not real*8 or integer') go to finish; endif; i=integers(ifbase,ifbase+isize-1); ibase=inewp; i=integers(ifbase,ifbase+isize-1); find=cc(i); call istrtor8(find,xx); x(j)=xx; if(inewp.eq.-99)go to finish; enddo; finish continue; return; end; == ==GLESJER Glesjer (1969) Heteroscedasticity Test subroutine glesjer(res,xdata,glesjert,prob,typetest,iprint1,iprint2); /; /; Performs the 15 variants of the Glesjer (1969) test for /; Heteroscedasticity. Ref /; Glesjer, H., "A New Test for Heteroscedasticity," /; Journal of the American Statistical Association, 64, 1969, pp 316-323 /; /; res = input Residual from equation to be tested /; xdata = input Columns to be used in the test. /; glesjer = 15 element array of Glesjer (1969) test statistics /; prob = 15 element array of probabilities of the Glesjer (1969) /; test. If prob(i) > .95 ther for that test honoscedasticity /; can be rejected. /; typetest= 15,2 names vector for the test to facitate Tables /; iprint1 = if set ne 0 prints internal regressions /; iprint2 = if set ne 0 prints test statistics /; Reference Greene (4th Edition page 510-511) /; /; Test involves three regresisons: /; /; 1. (res**2) = f(constant xdata) /; 2. |res| = f(constant xdata) /; 3. log(|res|) = f(constant xdata) /; /; A Wald test is used to test the joint restriction that all /; coefficients of the xdata are 0.0. The B34S Wald subroutine /; uses an F test. If this test is squared, it becomes a Chisquare /; test with two degreees of Freedom. /; /; Greene notes on page 511 "Since each of these regressions is /; heteroscedastic, we would use White's estimator of the covariance /; matrix in each case." However his test problem uses only the OLS /; estimate of the covariance. This implementation produces OLS, White, /; White1, White2, White3 results for each of the 3 formulations. /; /; Example: /; /; b34sexec options ginclude('greene4.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; /; /; xdata=catcol(%x(,3),%x(,4)); /; call glesjer(%res,xdata,glesjer,prob,typetest,iprint1, /; iprint2); /; /; b34srun; /; /; /; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /; glesjert=array(15:); prob =array(15:); if(iprint1.eq.0)print1=c1array(:' '); if(iprint1.ne.0)print1=c1array(:':print '); iprint=0; if(iprint2.ne.0)iprint=1; typetest=c8array(15,2:); t1=c8array(5:'OLS ','White ','White_1 ','White_2','White_3'); t2=c8array(3:'RES**2 ','|RES| ','log|RES|'); i=integers(1:5); ii=i; do jj=1,3; typetest(ii,1)=t2(jj); ii=ii+5; enddo; i=integers(1,5); ii=i; do jj=1,3; typetest(ii,2)=t1(i); ii=ii+5; enddo; q=vector(nocols(xdata):); r=matrix(nocols(xdata),nocols(xdata):)+1.0; i=nocols(xdata)+1; r(,i)=q; abs_e =dabs(res); e_sq =afam(res)*afam(res); lg_abs_e =dlog(abs_e); call olsq(e_sq xdata argument(print1) :qr ); if(iprint2.ne.0) call fprint(:clear :col 1 :display ' ' :print :clear :col 1 :display 'Glesjer (1969) Heteroscedasticity Test using RES**2.' :col 58 :display 'Variance-Covariance Matrix OLS.' :print); call wald(%coef,%xpxinv, r,q,%nob,%resvar,waldtest,probwald,iprint); glesjert(1)=waldtest; prob(1) =probwald; call olsq(e_sq xdata :white argument(print1) :qr); if(iprint2.ne.0) call fprint(:clear :col 1 :display ' ' :print :clear :col 1 :display 'Glesjer (1969) Heteroscedasticity Test using RES**2.' :col 58 :display 'Variance-Covariance Matrix White.' :print); call wald(%coef,%varcov2,r,q,%nob,1.0 ,waldtest,probwald,iprint); glesjert(2)=waldtest; prob(2) =probwald; call olsq(e_sq xdata :white1 argument(print1) :qr); if(iprint2.ne.0) call fprint(:clear :col 1 :display ' ' :print :clear :col 1 :display 'Glesjer (1969) Heteroscedasticity Test using RES**2.' :col 58 :display 'Variance-Covariance Matrix White Variant 1.' :print); call wald(%coef,%varcov2,r,q,%nob,1.0 ,waldtest,probwald,iprint); glesjert(3)=waldtest; prob(3) =probwald; call olsq(e_sq xdata :white2 argument(print1) :qr); if(iprint2.ne.0) call fprint(:clear :col 1 :display ' ' :print :clear :col 1 :display 'Glesjer (1969) Heteroscedasticity Test using RES**2.' :col 58 :display 'Variance-Covariance Matrix White Variant 2.' :print); call wald(%coef,%varcov2,r,q,%nob,1.0 ,waldtest,probwald,iprint); glesjert(4)=waldtest; prob(4) =probwald; call olsq(e_sq xdata :white3 argument(print1) :qr) ; if(iprint2.ne.0) call fprint(:clear :col 1 :display ' ' :print :clear :col 1 :display 'Glesjer (1969) Heteroscedasticity Test using RES**2.' :col 58 :display 'Variance-Covariance Matrix Variant 3.' :print); call wald(%coef,%varcov2,r,q,%nob,1.0 ,waldtest,probwald,iprint); glesjert(5)=waldtest; prob(5) =probwald; call olsq(abs_e xdata argument(print1) :qr); if(iprint2.ne.0) call fprint(:clear :col 1 :display ' ' :print :clear :col 1 :display 'Glesjer (1969) Heteroscedasticity Test using |RES|.' :col 58 :display 'Variance-Covariance Matrix OLS.' :print); call wald(%coef,%xpxinv, r,q,%nob,%resvar,waldtest,probwald,iprint); glesjert(6)=waldtest; prob(6) =probwald; call olsq(abs_e xdata :white argument(print1) :qr); if(iprint2.ne.0) call fprint(:clear :col 1 :display ' ' :print :clear :col 1 :display 'Glesjer (1969) Heteroscedasticity Test using |RES|.' :col 58 :display 'Variance-Covariance Matrix White.' :print); call wald(%coef,%varcov2,r,q,%nob,1.0 ,waldtest,probwald,iprint); glesjert(7)=waldtest; prob(7) =probwald; call olsq(abs_e xdata :white1 argument(print1) :qr); if(iprint2.ne.0) call fprint(:clear :col 1 :display ' ' :print :clear :col 1 :display 'Glesjer (1969) Heteroscedasticity Test using |RES|.' :col 58 :display 'Variance-Covariance Matrix White Variant 1.' :print); call wald(%coef,%varcov2,r,q,%nob,1.0 ,waldtest,probwald,iprint); glesjert(8)=waldtest; prob(8) =probwald; call olsq(abs_e xdata :white2 argument(print1) :qr); if(iprint2.ne.0) call fprint(:clear :col 1 :display ' ' :print :clear :col 1 :display 'Glesjer (1969) Heteroscedasticity Test using |RES|.' :col 58 :display 'Variance-Covariance Matrix White Variant 2.' :print); call wald(%coef,%varcov2,r,q,%nob,1.0 ,waldtest,probwald,iprint); glesjert(9)=waldtest; prob(9) =probwald; call olsq(abs_e xdata :white3 argument(print1) :qr) ; if(iprint2.ne.0) call fprint(:clear :col 1 :display ' ' :print :clear :col 1 :display 'Glesjer (1969) Heteroscedasticity Test using |RES|.' :col 58 :display 'Variance-Covariance Matrix Varient 3.' :print); call wald(%coef,%varcov2,r,q,%nob,1.0 ,waldtest,probwald,iprint); glesjert(10)=waldtest; prob(10) =probwald; call olsq(lg_abs_e xdata argument(print1) :qr); if(iprint2.ne.0) call fprint(:clear :col 1 :display ' ' :print :clear :col 1 :display 'Glesjer (1969) Heteroscedasticity Test using log(|RES|).' :col 58 :display 'Variance-Covariance Matrix OLS.' :print); call wald(%coef,%xpxinv, r,q,%nob,%resvar,waldtest,probwald,iprint); glesjert(11)=waldtest; prob(11) =probwald; call olsq(lg_abs_e xdata :white argument(print1) :qr); if(iprint2.ne.0) call fprint(:clear :col 1 :display ' ' :print :clear :col 1 :display 'Glesjer (1969) Heteroscedasticity Test using log(|RES|).' :col 58 :display 'Variance-Covariance Matrix White.' :print); call wald(%coef,%varcov2,r,q,%nob,1.0 ,waldtest,probwald,iprint); glesjert(12)=waldtest; prob(12) =probwald; call olsq(lg_abs_e xdata :white1 argument(print1) :qr); if(iprint2.ne.0) call fprint(:clear :col 1 :display ' ' :print :clear :col 1 :display 'Glesjer (1969) Heteroscedasticity Test using log(|RES|).' :col 58 :display 'Variance-Covariance Matrix White variant 1.' :print); call wald(%coef,%varcov2,r,q,%nob,1.0 ,waldtest,probwald,iprint); glesjert(13)=waldtest; prob(13) =probwald; call olsq(lg_abs_e xdata :white2 argument(print1) :qr); if(iprint2.ne.0) call fprint(:clear :col 1 :display ' ' :print :clear :col 1 :display 'Glesjer (1969) Heteroscedasticity Test using log(|RES|).' :col 58 :display 'Variance-Covariance Matrix White variant 2.' :print); call wald(%coef,%varcov2,r,q,%nob,1.0 ,waldtest,probwald,iprint); glesjert(14)=waldtest; prob(14) =probwald; call olsq(lg_abs_e xdata :white3 argument(print1) :qr) ; if(iprint2.ne.0) call fprint(:clear :col 1 :display ' ' :print :clear :col 1 :display 'Glesjer (1969) Heteroscedasticity Test using log(|RES|).' :col 58 :display 'Variance-Covariance Matrix White variant 3.' :print); call wald(%coef,%varcov2,r,q,%nob,1.0 ,waldtest,probwald,iprint); glesjert(15)=waldtest; prob(15) =probwald; return; end; == ==GLM_INFO Calculates out-of-sample best GLM model subroutine glm_info(y,holdout,glmcoef,glmconst,alm,glmf,fss,mod,iprint); /; /; inputs /; y => y variable from model /; holdout => holdout from model estimation /; glmcoef => %coef vector for the lmu models actuially run /; glmconst => %a0 vector of constants of the glm model /; lamda => %alm. Lamda value for model /; /; outputs /; glmf => glm forecasts out of sample for best model /; fss => forecast sum of squares /; mod => Model pointer /; iprint => =1 implies display a table /; /; Built 2 November 2009 /; n_to_test=nocols(glmcoef); fss=.1d+30; mod=1; do i=1,n_to_test; glmf=mfam(holdout)*vfam(glmcoef(,i)) + vfam(glmconst(i)); test=sumsq(vfam(glmf)-vfam(y)); if(iprint.ne.0) call fprint(:clear :string 'Model ' :col 7 :display i '(i4)' :col 12 :string 'RSS out of Sample ' :col 36 :display test '(g16.8)' :col 54 :string 'Lamda ' :col 60 :display alm(i) '(g16.8)' :print); if(test.lt.fss)then; mod=i; fss=test; endif; enddo; glmf=mfam(holdout)*vfam(glmcoef(,mod)) + vfam(glmconst(mod)); return; end; == ==GLMSOLVE Effect of Elastic Net Parameter on RSS & RSQ subroutine glmsolve(lhs,rhs,alower,aupper,na,nlam, thr,lam_min,thr,ne,rss,rsq,avalue,lam_val,holdout,iprint,igraph) /; /; Illustrate effect of Elastic Net Parameter on GLS solution /; /; lhs => Character Variable for left hand side. /; rhs => Character variable for right hand side /; use call makeglobal(all variables); /; call makelocat( ); /; alower => Lower Elastic Net Parameter /; aupper => Upper Elastic New Parameter /; na => Number of Elastic Net Parameters to calculate /; nlam => Number of Lamda values to consider /; thr => Convergence threshold /; ne => Sets maximum number of variables allowed in model /; rss => Residual sum of squares given avalue(i) /; rsq => R**2 given avalue(i) /; avalue => Calculated Elastic Net Parameter /; lam_val => Lamda value of best model given avalue /; holdout => Sample holdout. Can be 0 /; iprint => =0 no print /; =1 print tabulation only /; =2 print models estimated including OLS /; igraph => =0 no graphs produced /; =1 produce graphs rss_ep.wmf and rsq_ep.wmf /; /; Built November 2009 by Houston H. Stokes /; ------------------------------------------------------------- if(iprint.ne.0)call olsq(argument(lhs) argument(rhs):l1 :minimax :print); rss =array(na:); rsq =array(na:); avalue =array(na:); lam_val=array(na:); ajump=(aupper-alower)/dfloat(na-1); pp=array(40:)+1.; do i=1,na; avalue(i)=alower+(dfloat(i-1))*ajump; if(iprint.gt.1)call glm(argument(lhs) argument(rhs) :print :lamdamin lam_min :nlam nlam :holdout holdout /; :rel_penalty pp :thr thr :parm sfam(avalue(i)) :ne ne); if(iprint.le.1)call glm(argument(lhs) argument(rhs) :lamdamin lam_min :nlam nlam :holdout holdout /; :rel_penalty pp :thr thr :parm sfam(avalue(i)) :ne ne); k=imin(%rss); lam_val(i)=%alm(k); rss(i)=%rss(k); rsq(i)=%rsq(k); if(iprint.gt.1)call print(' ':); enddo; if(iprint.gt.0)call tabulate(avalue,rss,rsq,lam_val); if(igraph.ne.0)then; call graph(avalue,rss :plottype xyplot :nocontact :pgborder :file 'rss_ep.wmf' :grid :heading 'Residual sumsq vs Elastic Net Parameter'); call graph(avalue,rsq :plottype xyplot :nocontact :pgborder :file 'rsq_ep.wmf' :grid :heading 'R**2 vs Elastic Net Parameter'); endif; return; end; == ==GLS Cochrane & Orcutt GLS Estimation of AR(k) Model /; GLS does a Cochrane & Orcutt (1949) estimation of AR(k) Model /; Contains Program GLSMOD which determines the model /; Program GLSSET which sets defaults /; Program GLS which does the estimation. /; /; Note: GLS_ML estimated as ML AR(1) model using Rats Formula /; /; call glsset; => Placed right after call olsq( :savex); /; call gls; => Placed after call glsset; but after /; any custom switches /; /; Current glsset "defaults" /; /; %maxgls=1; /; %nl2sol=0; /; %plot=0; /; %maxit=2000; /; %flam = 1.0; /; %flu = 10.; /; %eps2 = .1e-12; /; /; /; Prior to call user sets /; %maxgls=k => Max GLS default = 1; /; %plot=1 => get residual plots /; %nl2sol=1 => Estimate with nl2sol in addition to NLLSQ /; /; Convergence switches /; /; %flam = 1.0; /; %flu = 10.; /; %eps2 = .1e-12; /; /; Example of a simple setup: /; 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=0; /; %plot=1; /; /; call gls; /; /; b34srun; /; /; ******************************************************************* /; Command built 1 April 2004 /; Estimation is GLS using nonlinear least squares. ML estimation for /; models is not possible for AR(k) models where k > 2. Nonlinear GLS /; Models can be estimated for any order in theory. /; /; This implementation drops the first observation. This is simular to /; the Cochrane and Orcutt approach (1949). For details see /; Greene (2003) page 273. Users wanting this adjustment can easily /; change the code. Most users will not require nl2sol which is provided /; as an option. Results obtained with this procedure can be compared to /; the REGRESSION command which does a two pass approach /; /; GLS test problem in staging.mac shows RATS AR1 test cases /; ******************************************************************* program glsmod; /$ /$ Form Yhat if GLS of order k /$ /$ Data in xp where each row is a series. /$ Parms in parm. %maxglas must be set /$ if(%maxgls.eq.0)then; yhat =parm*xp; endif; if(%maxgls.gt.0)then; p1=parm(integers(1,norows(parm)-%maxgls) ); p2=vector(:parm(integers( norows(parm)-%maxgls+1,norows(parm)))); yhat=p1*xp(,i); do j=1,%maxgls; yhat=yhat+(p2(j)* %y(i-j)); pp=p2(j)*p1; yhat=yhat-(pp*xp(,i-j)); enddo; endif; glsres=glsy-yhat; return; end; program gls; /; call echooff; %basex=%x; %basey=%y; if(kind(%maxgls).ne.-4)%maxgls=1; i=vfam(integers(%maxgls+1,norows(%y))); %olsres =%res; %olscoef=%coef; %olsse =%se; %olst =%t; if(%maxgls.gt.0)then; do kk=1,%maxgls; ii=norows(%olscoef)+kk; %names(ii)=object(rho_,eval('kk')); enddo; endif; parm=vector(norows(%coef)+%maxgls:)+.1; glsy=%y(i); xp=transpose(%x); /; save a copy to avoid "internal use of same name" maxit%=%maxit; flam%=%flam; flu%=%flu; eps2%=%eps2; call nllsq(glsy yhat :name glsmod :parms parm :maxit maxit% :flam flam% :flu flu% :eps2 eps2% :print result ); call tabulate(%names %lag %olscoef %olsse %olst %coef %se %t :title 'OLS and GLS Estimates using NLLSQ'); if(%plot.eq.1)call graph(%res :Heading 'GLS Residual NLLSQ' :nocontact :nolabel /; :markpoint 1 1 3 14 /; :colors black bblue :hardcopyfmt WMF :pspaceon :pgyscaleright 'i' :pgxscaletop 'i' :pgborder :file 'gls_res_nllsq.wmf' ); if(%nl2sol.eq.1)then; parm=vector(norows(%olscoef)+%maxgls:)+.1; call nl2sol(glsres :name glsmod :parms parm :maxit maxit% :ivalue parm :print /$ :itprint ); call tabulate(%names %lag %olscoef %olsse %olst %coef %se %t :title 'OLS and GLS Estimates using NL2SOL'); if(%plot.eq.1) call graph(%res :Heading 'GLS Residual using NL2SOL' :nocontact :nolabel /; :markpoint 1 1 3 14 /; :colors black bblue :hardcopyfmt WMF :pspaceon :pgyscaleright 'i' :pgxscaletop 'i' :pgborder :file 'gls_res_nl2sol.wmf' ); endif; return; end; program glsset; %maxgls=1; %nl2sol=0; %plot=0; %maxit=2000; %flam = 1.0; %flu = 10.; %eps2 = .1e-12; return; end; == ==GLSDATA Transform x to x* given GLS rho vector 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 /$ maxgls=norows(rho); i = integers(maxgls+1,norows(x)); ii = i-maxgls; if(nocols(x).gt.1)then; newx=x(i,); do j=1,maxgls; newx(ii,)=newx(ii,)-(sfam(rho(j))*x((i-j),)); enddo; endif; if(nocols(x).eq.1)then; newx=x(i); do j=1,maxgls; newx(ii)=newx(ii)-(sfam(rho(j))*x((i-j))); enddo; endif; return; end; == ==GLS_ML GLS_ML Maximum Likelihood estimation of AR(1) Model /; GLS_ML Maximum Likelihood estimation of AR(1) Model /; /; The alternative GLS does a Cochrane & Orcutt (1949) /; estimation of an AR(k) Model /; /; GLS_ML Contains Program GLSMLMOD which determines the model /; Program GLSMLSET which sets defaults /; Program GLS_ML which does the estimation. /; /; /; call glsmlset; => Placed right after call olsq( :savex); /; call gls_ml; => Placed after call glsmlset; but after /; any custom switches /; /; Current glsml_set "defaults" /; /; program glsmlset; /; %plot=0; /; _%maxf =2000; /; _%maxg =2000; /; _%maxit=2000; /; return; /; end; /; /; /; Prior to call user sets /; %plot=1 => get residual plots /; /; Example of a simple setup: /; b34sexec matrix; /; call loaddata; /; call load(gls_ML :staging); /; call olsq(gasout gasin{1 to nn} gasout{1 to nn} /; :print :savex); /; call glsmlset; /; %plot=1; /; /; call gls_mL; /; /; b34srun; /; /; ******************************************************************* /; Command built 1 May 2004 /; Estimation is GLS using ML. ML estimation for AR(k) models /; where k > 1 is not possible with this routine. /; /; Nonlinear GLS Models can in theory be estimated for any order. /; /; ML estimation for AR(k) models where k > 2 is not possible. /; /; Greene (2003, 274) notes that Beach and MacKinnon (1978) have /; derived an algorithm for AR(2) disturbances. Greene notes /; "For higher-order autoregressive models maximum likelihood /; estimation is presently impractical, but the two-step estimators /; can easily be extended." /; /; Nonlinear GLS Models can in theory be estimated for any order. /; /; The GLS_ML test problem in staging.mac shows RATS AR1 test cases /; The GLS_ML test problem in staging.mac shows RATS AR1 test cases /; ******************************************************************* program glsmlmod; /$ /$ Form function for use with ML AR(1) estimation /$ p1=parm(integers(1,norows(parm)-2)); p2=parm(norows(parm)-1); s2=parm(norows(parm)); yhat=p1*xp(,i); yhat=yhat+(p2* %y(i-1)); pp=p2*p1; yhat=yhat-(pp*xp(,i-1)); glsres=glsy-yhat; part1=.5*dlog(dmax1((1.0-(p2*p2)),.1e-8)); part2=((-.5)*(1.0-(p2*p2))*((%y(1)-xp(,1)*p1)**2.))/s2; part3= (-.5)*dfloat(norows(yhat))*dlog(dmax1(s2,.1e-8)); func =part1+part2+part3-(.5*(sumsq(glsres))/s2); return; end; program glsmlset; %plot=0; _%maxf =2000; _%maxg =2000; _%maxit=2000; return; end; program gls_ml; call echooff; i=vfam(integers(2,norows(%y))); %olsres =%res; %olscoef=%coef; %olsse =%se; %olst =%t; ii=norows(%olscoef)+1; %names(ii)=object(rho_,1); parm=vector(norows(%coef)+1:)+.1; jj=norows(parm)-1; rvec=parm; rvec(jj)=%coef(jj); call print(glsmlmod); s2=%resvar; glsy=%y(i); xp=transpose(%x); /; save a copy to avoid "internal use of same name" nparms=norows(parm); ll=array(nparms+1:)-100.0d+10; uu=array(nparms+1:)+100.0d+10; ll(nparms)=-.999; uu(nparms)= .999; npp=nparms+1; ll(npp)=.1e-14; parm(npp)=%resvar; rvec(npp)=%resvar; %names(npp)=object(s_,2); call cmaxf2(func :name glsmlmod :parms parm :ivalue rvec :lower ll :maxit _%maxit /$ :steptol .1e-15 :maxfun _%maxf :maxg _%maxg :upper uu :print); call tabulate(%names %lag %olscoef %olsse %olst %coef %se %t :title 'OLS and ML_GLS Estimates using CMAXF2'); %res=glsres; if(%plot.eq.1)then; call graph(%res :Heading 'ML GLS AR(1) Residual using CMAX2' :nocontact :nolabel /; :markpoint 1 1 3 14 /; :colors black bblue :hardcopyfmt WMF :pspaceon :pgyscaleright 'i' :pgxscaletop 'i' :pgborder :file 'gls_res_ml_cmax2.wmf' ); endif; return; end; == ==GLS_2P Two Pass Generalized Least Squares for Order K program gls2pset; /; Internal Default Settings %print =1; %psteps=0; %maxgls=1; %doplots=0; %doacf=0; return; end; program gls_2P; /; Calculated GLS Using Two Variants of Two Pass Method /; /; Method # A traps the residual vector from OLS and fits /; e(t)=rho(1)*e(t-1)+rho(2)*e(t-2)+...+rho(k)*e(t-k) /; /; Data Smoothed. and adjusted model estimated /; /; Method # B uses the Classic B34T/B34S approach that attempts to /; estimate each rho with the max amount of data /; /; /; Use: /; call load(glsdata :staging); /; call load(gls_2p) :staging); /; call olsq y x1 x2 :savex :print); /; call gls2pset; /; call gls_2p; /; /; Optional Switches placed after call gls2pset; /; %maxgls=1 ... Default=1 /; %print=1 or 0 Default = 1 /; %psteps =0 or 1; Default = 0 /; %doplots=1 Default = 0 => Plots raw and GLS res /; %doacf =0 Default = 0 => # of ACF to calculate /; Must be set if plots requested /; /; Settings in gls2pset /; /; program gls2pset; /; %print =1; /; %psteps=0; /; %maxgls=1; /; %doacf=0; /; %doplots=0; /; return; /; end; /; /; Variables Created: /; /; %coef_a => GLS Coef Method A /; %se_a => GLS SE Method A /; %t_a => GLS T Method A /; %coef_b => GLS Coef Method B /; %se_b => GLS SE Method B /; %t_b => GLS T Method B /; %coef => OLS Coef /; %se => OLS SE /; %t => OLS T /; %baseres => OLS Residual /; %basey => OLS Y varuable /; %res_a => GLS Residuale Method A /; %y_a => GLS Method A Y variable /; %x_a => GLS Method A X variables /; %res_b => GLS Residual Method B /; %y_b => GLS Method B Y variable /; %x_b => GLS Method A X variables /; %rho_a => Rho for last step using method A /; %rho_b => Rho for last step using method B /; /; call free(%acfols,%acfres_a,%acfres_b); %basenob=%nob; %basex=%x; %basey=%y; %baseres=%res; %holdnam=%names; %c_base=%coef; %sebase=%se; %tbase =%t; xpx=matrix(%maxgls,%maxgls:); xpy=vector(%maxgls:); do jj=1,%maxgls; %print2=0; if(%psteps.eq.1 .and.%print.eq.1)%print2=1; if(jj.eq.%maxgls.and.%print.eq.1)%print2=1; if(%print2.eq.1) call olsq(%baseres %baseres{1 to jj} :noint :print :savex); if(%print2.ne.1) call olsq(%baseres %baseres{1 to jj} :noint :savex); %xhold=%x; %yhold=%y; call glsdata(%coef,%basey,%y_a); call glsdata(%coef,%basex,%x_a); %rho_a=%coef; if(%print2.eq.1)then; call olsq(%y_a %x_a :noint :print); call print(' ':); call print('+++++++++++++++++++++++++++++++++++++++++++++++++':); call print('GLS Rho Coefficients for Variant A',%rho_a); call print('+++++++++++++++++++++++++++++++++++++++++++++++++':); call print(' ':); endif; if(%print2.ne.1)call olsq(%y_a %x_a :noint ); %coef_a=%coef; %se_a =%se; %t_a =%t; %res_a=%res; xpx2=transpose(%xhold)*%xhold; xpy2=transpose(%xhold)*%yhold; do kk=1,jj; do mm=jj,jj; xpx(kk,mm)=xpx2(kk,mm); if(kk.ne.mm)xpx(mm,kk)=xpx2(kk,mm); enddo; enddo; xpy(jj)=xpy2(jj); /; Define sum matrix; xpxnew=submatrix(xpx,1,jj,1,jj); xpynew=xpy(integers(1,jj)); /; adjust vectors for DF if(klass(xpynew).eq.0)xpynew=vector(:xpynew); do iadj=1,jj; do jadj= iadj,jj; xpxnew(iadj,jadj)= xpxnew(iadj,jadj)*(1.0/dfloat(%basenob-jadj)); if(iadj.ne.jadj)xpxnew(jadj,iadj)=xpxnew(iadj,jadj); enddo; xpynew(iadj)= xpynew(iadj)* (1.0/dfloat(%basenob-iadj)); enddo; %rho_b=inv(xpxnew)*xpynew; if(%print2.eq.1)then; call print('GLS Rho Coefficients using Variant # B',%rho_b); call print('+++++++++++++++++++++++++++++++++++++++++++++++++':); endif; call glsdata(%rho_b,%basey,%y_b); call glsdata(%rho_b,%basex,%x_b); if(%print2.ne.1)call olsq(%y_b %x_b :noint ); if(%print2.eq.1)call olsq(%y_b %x_b :noint :print); %coef_b= %coef; %se_b = %se; %t_b = %t; %res_b = %res; %coef = %c_base; %se = %sebase; %t = %tbase; if(%print2.eq.1)then; call print(' ':); call print('GLS Using Two Pass Method Variants A & B':); call tabulate(%holdnam %coef %se %t %coef_a %se_a %t_a %coef_b %se_b %t_b); endif; enddo; if(%doplots.eq.1.and.%doacf.ne.0)then; call graph(%res_a,%res_b); call character(cc,'Residuals from OLS Model'); call data_acf(%baseres,cc,%doacf); call character(cc,'Residuals from GLS Model using Method A'); call data_acf(%res_a,cc, %doacf); call character(cc,'Residuals from GLS Model using Method B'); call data_acf(%res_b,cc, %doacf); endif; if(%doacf.ne.0)then; %acfols =acf(%baseres,%doacf); %acfresa=acf(%res_a, %doacf); %acfresb=acf(%res_b, %doacf); call tabulate(%acfols,%acfresa,%acfresb); if(%doplots.ne.0)call graph(%acfols,%acfresa,%acfresb :nolabel); endif; return; end; == ==GPH Geweke & Porter-Hudak (1983) Fractional Diff. Test subroutine gph(series,power,d,se,se2,iprint); /; /; Estimates fractional difference power for series using the /; frequency domain regression techniques of Geweke and Porter-Hudak. /; /; Subroutine GPH works in a manner similar to the RATS GPH routine /; which was used as a basis to compare against. /; /; series => series to analyze. This is usually the first difference /; of the series actually being studied. /; /; POWER => Lowest T**power frequencies are used in running the /; regression. Usually set as .5 (i.e. use square root of the /; number of observations /; /; d => estimate of fractional differencing /; /; se => Standard error of d /; /; se2 => Asymptotic SE of d /; /; iprint => 0 => no printing /; 1 => print results /; 2 => print results and regression /; /; Reference: Geweke and Porter-Hudak, "The Estimation and Application /; of Long Memory Time Series Models", J. of Time Series Analysis, /; 1983, (221-238). /; /; Built 30 April 2006 by Houston H. Stokes. /; Improved error messages 5 September 2006. /; call spectral(series,sinx,cosx,px,sx,freq1); n=idint(dfloat(norows(series))**power); if(n.lt.4)then; call epprint('ERROR: Power too small in call gph( ). # obs LT 4'); call epprint(' (norows(series)**power) .lt. 4':); call stop; endif; if(n.gt.((norows(series)/2)-5))then; call epprint('ERROR: Power too large in call gph( ). '); call epprint(' (norows(series)**power)>(norows(series)/2)-5)':); call stop; endif; k=grid(1., dfloat(n), 1.); ftest=2. * pi() * k / dfloat(norows(series)); t_freq=log(4. * (sin(afam(freq1)/2.)**2.)); y = log(px(integers(1,n))); x = t_freq(integers(1,n)); if(iprint.eq.2)call olsq(y x :print); if(iprint.ne.2)call olsq(y x ); se=%se(1); se2= pi()*sqrt(%xpxinv(1,1)/6.0); d=(-1.)*%coef(1); z =d/se; z2=d/se2; if(iprint.ne.0)then; call print(' ':); call print('Geweke & Porter-Hudak (1983) Fractional Difference Test':); call print('Power ',power:); call print('Number of Frequencies used ',%nob :); call print('Estimated d ',d :); call print('Asymptotic Standard Error ',se2 :); call print('OLS Standard Error ',se :); call print('Asymptotic Z score ',z2 :); call print('OLS t score ',z :); endif; return; end; == ==G_QUANDT Goldfeld-Quandt Heteroscedasticity Test subroutine g_quandt(y,x,sortcol,name,lag1,n1,n2,ss1,ss2,g_q_test,prob, iprint1,iprint2); /; /; 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 /; name => Names Col /; lag1 => 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 /; ++++++++++++++++++++++++++++++++++++++++++++++++++++ /; /; Example from Greene (4th Ed Page 510): /; b34sexec options ginclude('greene4.mac') member(a5_1); /; b34srun; /; /; b34sexec matrix; /; call loaddata; /; call load(g_quandt :staging); /; call echooff; /; mmm =(exp .gt. 0.0); /; call olsq(exp age ownrent income incomesq /; :sample mmm :print :savex); /; /; do i=3,3; /; call g_quandt(%y,%x,i,%names,%lag,36,37,ss1,ss2,g_q,prob,1,1); /; enddo; /; b34srun; /; /; +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /; Command built December 2005 /; o_order=integers(1,norows(x)); yuse=y; xuse=x; if(sortcol.ne.0)then; ii=ranker(x(,sortcol)); yuse=y(ii); xuse=x(ii,); endif; m1=(o_order.le.n1); m2=(o_order.ge.n2); if(iprint1.ne.0)then; call olsq(yuse xuse :noint :sample m1 :print :qr); ss1 =%rss; call olsq(yuse xuse :noint :sample m2 :print :qr); ss2 =%rss; endif; if(iprint1.eq.0)then; call olsq(yuse xuse :noint :sample m1 :qr); ss1 =%rss; call olsq(yuse xuse :noint :sample m2 :qr); ss2 =%rss; endif; g_q_test=ss1/ss2; df1 =dfloat(n1-%k); df2 =dfloat(norows(y)-n1-%k); g_q_test=g_q_test*(df2/df1); if(g_q_test.ge.1.0d+00)prob=fprob(g_q_test, df1,df2); if(g_q_test.lt.1.0d+00)then; g_q_test=1.0/g_q_test; prob=fprob(g_q_test,df2,df1); endif; if(iprint2.ne.0)then; call print(' ':); if(sortcol .ne.0)then; work=catrow( 'Data sorted by ', c1array(8:name(sortcol))); call print(c1array(:work):); if(lag1(sortcol) .ne.0) call print('Lag was ', lag1(sortcol) ); endif; isdrop=n2-n1-1; call fprint(:clear :display ' ' :print :col 1 :display 'Sample 1 ends at obs ' :col 30 :display n1 '(i5)' :print :clear :col 1 :display 'Sample 2 starts from obs' :col 30 :display n2 '(i5)' :print :clear :col 1 :display '# dropped ' :col 30 :display isdrop '(i5)' :print :clear :col 1 :display 'Goldfeld-Quandt (1965) Heteroscedasticity Test ' :col 48 :display g_q_test '(g16.8)' :col 65 :display 'Probability' :col 77 :display prob '(g12.4)' :print); endif; return; end; == ==HETTEST1 Theil Heteroskedasticity Test on Sorted Data subroutine hettest1(res,x,name,lags,f,rf,sig, isort1,isort2,pertest,print1); /; /; res => Residual vector usually %res /; x => x matrix %x /; name => Used for Goldfeld-Quandt test /; lags => Used for Goldfeld-Quandt test /; f => first pertest squared / last pertest squared /; rf => 1/f /; sig => Significance of f or rf /; isort1 => Beginning col to sort for Het Test /; isort2 => End col to sort for Het Test /; if isort1 or isort2 = 0 => no sorting is done. /; pertest => set in range .1 - .5 => sets # dropped /; print1 => print test /; /; Routine built January 2008. Arguments subject to change /; Test based on Theil 1971 test that was developed for BLUS /; residuals /; /; /; /; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /; Example: /; /; 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(het_test1 :staging); /; /; mmm =(exp .gt. 0.0); /; call olsq(exp age ownrent income incomesq /; :sample mmm :print :savex); /; /; /; Parameters for het_test /; /; isort1=3; /; isort2=3; /; pertest=.33333; /; print1=1; /; /; call hettest1(%res,%y,%x,%names,%lag,f,rf,sig,isort1,isort2, /; pertest,print1); /; /; b34srun; /; /; ++++++++++++++++++++++++++++++ end of extended example ++++++++++ /; n=norows(x); /; /; p=pertest; if(p.gt..5)p=.5; n1=idnint(p*dfloat(n)); /; n1 = the number in top res(1) ... res(n1) /; n2 = bottom res(n2) ... res(n) /; nn2 = number in bottom n2=n-n1; nn2=n-n2+1; if((n1+nn2).gt.n)then; n2=n2+1; nn2=nn2-1; endif; /; if(isort1.eq.0.and.isort2.eq.0)then; top=sumsq(res(integers(1,n1))); bot=sumsq(res(integers(n2,n))); f =missing(); rf=missing(); if(bot.ne.0.0)f=top/bot; if(top.ne.0.0)rf=1./f; if(f.ge.1.0 )sig=fprob(f,dfloat(n1), dfloat(nn2)); if(rf.ge.1.0)sig=fprob(rf,dfloat(nn2),dfloat(n1)); if(print1.ne.0)then; call fprint(:col 1 :string 'F(' :col 3 :display n1 '(i6)' :col 9 :string ',' :col 10 :display nn2 '(i6)' :col 16 :string ')' :col 20 :display f '(g16.8)' :col 40 :string '1/F' :col 45 :display rf '(g16.8)' :col 62 :string 'Hetroskedasticity' :col 82 :display sig '(f8.4)' :col 92 :string 'Original Order of Data' :print :clear); endif; go to done; endif; if((isort1.gt.0.and.isort2.le.0).or.(isort1.le.0.and.isort2.gt.0))then; call epprint('ERROR: either isort1 or isort2 =0':); call stop; endif; if(isort1.ne.0.and.isort2.ne.0)then; ntests=isort2-isort1+1; if(ntests.le.0)then; call epprint('ERROR: isort2 and isort1 imply ntests LE 0':); call stop; endif; f =array(ntests:); rf =array(ntests:); sig =array(ntests:); j=0; do i=isort1,isort2; if(i.gt.nocols(x))then; call epprint('ERROR: isort1 ands isort2 imply col gt that in x':); call stop; endif; j=j+1; localres=res(ranker(x(,i))); top=sumsq(localres(integers(1,n1))); bot=sumsq(localres(integers(n2,n))); f(j) =missing(); rf(j)=missing(); if(bot.ne.0.0)f(j)=top/bot; if(top.ne.0.0)rf(j)=1./f(j); if(f(j).ge. 1.0)sig(j)=fprob(f(j), dfloat(n1), dfloat(nn2)); if(rf(j).ge.1.0)sig(j)=fprob(rf(j),dfloat(nn2),dfloat(n1)); if(print1.ne.0)then; call fprint(:col 1 :string 'F(' :col 3 :display n1 '(i6)' :col 9 :string ',' :col 10 :display nn2 '(i6)' :col 16 :string ')' :col 20 :display f(j) '(g16.8)' :col 40 :string '1/F' :col 45 :display rf(j) '(g16.8)' :col 62 :string 'Hetroskedasticity' :col 82 :display sig(j) '(f8.4)' :col 92 :string 'for' :col 96 :display name(i) :col 104 :string '{' :col 105 :display lags(i) '(i4)' :col 109 :string '}' :print :clear); endif; enddo; done continue; return; end; == ==HET_TEST Heteroscedasticity Testing subroutine het_test(res,y,x,name,lags,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 /; name => 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 /; gqprint => 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 /; /; Requires: call load(b_g_test); /; /; /; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /; Example: /; /; 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; /; /; /; Parameters for het_test /; /; n1=36; /; n2=37; /; gqprint1=0; /; gqprint2=1; /; isort1=3; /; isort2=3; /; /; 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; /; /; ++++++++++++++++++++++++++++++ end of extended example ++++++++++ /; yhat=y-res; /; Breusch-Pagan Test Three variants call print(' ':); holdres=res; holdnob=norows(res); holdx=x; sig=sumsq(holdres)/dfloat(holdnob); adj_res=(afam(holdres)*afam(holdres))/sig; z=dfloat(integers(1,holdnob)); /; /; LM test appears to be same as against t /; /; /; g=vfam(adj_res-1.0); /; lm=.5*g*x*inv(transpose(x)*x)*transpose(x)*g; /; df=2.0; /; call fprint(:clear :display 'Breusch-Pagan LM Het. test ' /; :col 35 :display lm '(g12.4)' /; :col 50 :display 'Chi-Sqr' /; :col 57 :display df '(f6.0)' /; :col 65 :display 'Probability' /; :col 77 :display chisqprob(lm,df) '(g12.4)' /; :print); if(bp1.ne.0)then; call olsq(adj_res z :qr); test=(((variance(adj_res)*dfloat(norows(adj_res)-1)))-%rss)/2.0; df=1.0; call fprint(:clear :display 'Breusch-Pagan Het. test using t ' :col 35 :display test '(g12.4)' :col 50 :display 'Chi-Sqr' :col 57 :display df '(f6.0)' :col 65 :display 'Probability' :col 77 :display chisqprob(test,1.0) '(g12.4)' :print); endif; if(bp2.ne.0)then; call olsq(adj_res holdx :noint :qr); test=(((variance(adj_res)*dfloat(norows(adj_res)-1)))-%rss)/2.0; df=dfloat(nocols(holdx)-1); call fprint(:clear :display 'Breusch-Pagan Het. test using X ' :col 35 :display test '(g12.4)' :col 50 :display 'Chi-Sqr' :col 57 :display df '(f6.0)' :col 65 :display 'Probability' :col 77 :display chisqprob(test,df) '(g12.4)' :print); endif; newx=array(holdnob,nocols(holdx)*2-1:); do i=1,nocols(holdx); newx(,i)=afam(holdx(,i)); enddo; j=1; do i=nocols(holdx)+1,nocols(newx); newx(,i)=(afam(holdx(,j))*afam(holdx(,j))); j=j+1; enddo; /; call print(holdx,newx); if(bp3.ne.0)then; call olsq(adj_res newx :noint :qr); test=(((variance(adj_res)*dfloat(norows(adj_res)-1)))-%rss)/2.0; df=dfloat(nocols(newx)-1); call fprint(:clear :display 'Breusch-Pagan Het. test using X**2' :col 35 :display test '(g12.4)' :col 50 :display 'Chi-Sqr' :col 57 :display df '(f6.0)' :col 65 :display 'Probability' :col 77 :display chisqprob(test,df) '(g12.4)' :print); endif; /; White (1980) Het Test if(white1.ne.0)then; call print(' ':); adj_res=(afam(holdres)*afam(holdres)); z=dfloat(integers(1,holdnob)); call olsq(adj_res z :qr); test=dfloat(%nob)*%rsq; df=1.0; call fprint(:clear :display 'White (1980) Het. test using t ' :col 35 :display test '(g12.4)' :col 50 :display 'Chi-Sqr' :col 57 :display df '(f6.0)' :col 65 :display 'Probability' :col 77 :display chisqprob(test,1.0) '(g12.4)' :print); endif; if(white2.ne.0)then; call olsq(adj_res holdx :noint :qr); test=dfloat(%nob)*%rsq; df=dfloat(nocols(holdx)-1); call fprint(:clear :display 'White (1980) Het. test using X ' :col 35 :display test '(g12.4)' :col 50 :display 'Chi-Sqr' :col 57 :display df '(f6.0)' :col 65 :display 'Probability' :col 77 :display chisqprob(test,df) '(g12.4)' :print); call fprint(:clear :display 'Same as Rats version of Koenker-Bassett (1982) Breusch-Pagan Test' :print); endif; if(white3.ne.0)then; newx=array(holdnob,nocols(holdx)*2-1:); do i=1,nocols(holdx); newx(,i)=afam(holdx(,i)); enddo; j=1; do i=nocols(holdx)+1,nocols(newx); newx(,i)=(afam(holdx(,j))*afam(holdx(,j))); j=j+1; enddo; call olsq(adj_res newx :noint :qr); test=dfloat(%nob)*%rsq; df=dfloat(nocols(newx)-1); call fprint(:clear :display 'White (1980) Het. test using X**2' :col 35 :display test '(g12.4)' :col 50 :display 'Chi-Sqr' :col 57 :display df '(f6.0)' :col 65 :display 'Probability' :col 77 :display chisqprob(test,df) '(g12.4)' :print); endif; /; +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ if(n1.ne.0.and.n2.ne.0)then; if(isort1.ne.0.and.isort2.ne.0)then; do i=isort1,isort2; call g_quandt(y,x,i,name,lags,n1,n2,ss1,ss2,g_q,prob, gqprint1,gqprint2); enddo; endif; if(isort1.eq.0.and.isort2.eq.0) call g_quandt(y,x,0,name,lags,n1,n2,ss1,ss2,g_q,prob, gqprint1,gqprint2); endif; return; end; == ==IACF Inverse ACF subroutine iacf(series,iacf1,n); /; /; Inverse ACF Experimental !!!!!!!!! /; /; Assume ARMA(p,q) model. Estimate ACF of ARMA(q,p) model /; /; series => series to estimate /; iacf1 => Inverse ACF /; n => N of terms. If n set too small /; AR(n) will not filter series /; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /; /; Routine built 7 June 2006 /; Routine uses logic from Wei /; x=series-mean(series); call olsq(x x{1 to n} :noint :qr); base=sumsq(%coef)+1.0; iacf1=array(n:); do i=1,n; top=-1.*%coef(i); do j=1,n-i+1; if((j+i).le.n)then; top=top+%coef(j)*%coef(j+i); endif; enddo; iacf1(i)=top/base; enddo; return; end; == ==IACFN Inverse Autocorrelation subroutine iacfn(series,iacf1,n); /; /; Inverse ACF Experimental - uses pacf following /; 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 /; /; Assume ARMA(p,q) model. Estimate ACF of ARMA(q,p) model /; /; series => series to estimate /; iacf1 => Inverse ACF /; n => N of terms. If n set too small /; AR(n) will not filter series /; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /; /; Routine built 7 June 2006 /; Bug fixed 16 July 2007 /; x=series-mean(series); acf1=acf(x,n,se,pacf); base=sumsq(pacf)+1.; iacf1=array(n:); do i=1,n; top=-1.*pacf(i); do j=1,n-i+1; if((j+i).le.n)then; top=top+(pacf(j)*pacf(j+i)); endif; enddo; iacf1(i)=top/base; enddo; return; end; == ==LASSO LASSO Shrinkage Approach subroutine lasso(y,x2,olscoef,lcoef,l_t,lamda,lresid,iprint); /; /; For discussion see Stokes (200x, chapter 10) /; /; Implements the LASSO shrinkage Method /; Reference: Hastie-Tibshirani-Friedman (2001) Page 64 and 77 /; See also page 72 for generalized setup which was used in /; this implimentation /; /; y => left hand side /; x => Right hand side. Usually %x from :savex /; olscoef => %coef from call olsq. Used for starting values /; lcoef => Lasso Coef /; l_t => lasso Coef t /; lamda => Lamda for Lasso Model. Larger Lamda => more shrinkage /; lresid => Residual from Lasso Model /; iprint => =0 No print use cmax2 /; iprint => =0 Print use cmax2 /; iprint => =0 No print use max2 /; iprint => =0 Print use max2 /; /; Note: The constant is not restricted!! /; Page 77 illustrates a centered example /; +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rvec=array(:olscoef); ll= array(norows(rvec):)-1.d+32; uu= array(norows(rvec):)+1.d+32; x=mfam(x2); lcoef=vfam(olscoef); i=integers(norows(olscoef)-1); if(iprint.eq.0) call cmaxf2(func :name lasso_1 :parms lcoef :ivalue rvec :maxit 2000 :lower ll :upper uu); if(iprint.eq.1) call cmaxf2(func :name lasso_1 :parms lcoef :ivalue rvec :maxit 2000 :lower ll :upper uu :print); if(iprint.eq.2) call maxf2(func :name lasso_1 :parms lcoef :ivalue rvec :maxit 2000); if(iprint.eq.3) call maxf2(func :name lasso_1 :parms lcoef :ivalue rvec :maxit 2000 :print); lresid =sumsq(y-x*lcoef); l_t=%t; lse=%se; if(iprint.eq.1.or.iprint.eq.3)then; call print('Lamda for Lasso model ',lamda:); call print('Sum of squared Residuals for Lasso Model ',lresid); endif; return; end; program lasso_1; func=(-1.0)*(sumsq(y-x*lcoef) + lamda*(sum(abs(lcoef(i)))) ); call outstring(3,3,'Function'); call outdouble(36,3,func); return; end; == ==LASSO2 Generalized LASSO Shrinkage Approach subroutine lasso2(y,x2,olscoef,lcoef,l_t,lamda,lresid,iprint,lamda_d); /; /; For discussion see Stokes (200x, chapter 10) /; /; Implements the LASSO shrinkage Method /; Reference: Hastie-Tibshirani-Friedman (2001) Page 64 and 77 /; See also page 72 for generalized setup which was used in /; this implimentation /; /; y => left hand side /; x => Right hand side. Usually %x from :savex /; olscoef => %coef from call olsq. Used for starting values /; lcoef => Lasso Coef /; l_t => lasso Coef t /; lamda => Lamda for Lasso Model. Larger Lamda => more shrinkage /; lresid => Residual from Lasso Model /; iprint => =0 No print use cmax2 /; iprint => =0 Print use cmax2 /; iprint => =0 No print use max2 /; iprint => =0 Print use max2 /; lamda_d => exponent on | | term /; =1.0 => lasso /; =2.0 => Ridge on Centered Data /; =0.0 => penalized subset OLS /; /; Note: The constant is not restricted!! /; Page 77 illustrates a centered example /; +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rvec=array(:olscoef); ll= array(norows(rvec):)-1.d+32; uu= array(norows(rvec):)+1.d+32; x=mfam(x2); lcoef=vfam(olscoef); i=integers(norows(olscoef)-1); if(iprint.eq.0) call cmaxf2(func :name lasso_2 :parms lcoef :ivalue rvec :maxit 2000 :lower ll :upper uu); if(iprint.eq.1) call cmaxf2(func :name lasso_2 :parms lcoef :ivalue rvec :maxit 2000 :lower ll :upper uu :print); if(iprint.eq.2) call maxf2(func :name lasso_2 :parms lcoef :ivalue rvec :maxit 2000); if(iprint.eq.3) call maxf2(func :name lasso_2 :parms lcoef :ivalue rvec :maxit 2000 :print); lresid =sumsq(y-x*lcoef); l_t=%t; lse=%se; if(iprint.eq.1.or.iprint.eq.3)then; call print('Lamda for Lasso model ',lamda:); call print('Lamda exponent ',lamda_d:); call print('Sum of squared Residuals for Generalized Lasso Model ', lresid); endif; return; end; program lasso_2; func=(-1.0)*(sumsq(y-x*lcoef) + lamda*(sum(afam(abs(lcoef(i)))**lamda_d))); call outstring(3,3,'Function'); call outdouble(36,3,func); return; end; == ==LTS Least Trimmed Squares subroutine lts(y,x,resid,names,lags,oldcoef,oldse,oldt,p,sort_y,sort_x, ihold,iprint); /; /; Least Trimmed Squares /; /; Reference: "Linear Models with R" By Julian Faraway (2005) /; Page 101 /; /; y => left hand side. Usually %y /; x => right hand side. Usually %x /; resid => Residual from original regression Usually %res /; names => Names from regression. Usually %names /; lags => lags from original regression. Usually %lags /; oldceof => Coefficient for full model. Usually %coef /; oldse => SE from original Model. Usually %se /; oldt => t from original model. Usually %t /; p => 1 - % trimmed /; sort_y => Y after sorting and truncation /; sort_x => X after sorting and truncation /; ihold => # of obs held out given p. /; iprint => If > 0 => print. /; /; /; Help on how to recover coef outside of routine /; Assume: /; /; 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); /; option # 1 /; call olsq(newy newx :noint :print :holdout ihold); /; option # 2 /; n=norows(newy); /; call deleterow(newy,n-ihold+1,ihold); /; call deleterow(newx,n-ihold+1,ihold); /; call print(inv(transpose(newx)*newx)*transpose(newx)*newy); /; /; ------------------------------------------------------------------- /; /; LTS is an example of a resistent regression method. The objective is /; to see how sensitive the results are to outliers. /; /; Built 21 June 2007 by Houston H. Stokes /; Mods 23 June 2007 /; n=norows(x); if(p.lt.0. .or. p .gt.1.)then; call epprint('p not in range 0 lt p le 1. Was ',p:); go to done; endif; r=afam(resid)*afam(resid); j=ranker(r); sort_x=x(j,); sort_y=y(j); ihold=idint((1.0-p)*dfloat(n)); if(iprint.eq.0) call olsq(sort_y sort_x :noint :qr :holdout ihold); if(iprint.ne.0)then; call print(' ':); call print('Least Trimmed Squares with holdout ',ihold:); call olsq(sort_y sort_x :noint :print :qr :holdout ihold); call print(' ':); call print('Least Trimmed Squares with holdout % ',(1.0-p):); call tabulate(names lags oldcoef oldse oldt %coef %se %t :cname); endif; done continue; return; end; == ==LTS_REC Recursive Least Trimmed Squares - Program program lts_rec; /; /; Does Recursive lts. Subroutine lts must be called first /; /; Needs the following set: /; /; oldcoef =1; /; n_recur =6; /; /; oldcoef = 0 => use Base Coefficients for table /; oldcoef = 1 => use prior LTS coef for table /; n_recur = => Sets number of recursive LTS estimates /; /; Warning: 1. Call lts must be called with exactly these argument /; names. /; /; 2. newy, newx and ihold are changed /; /; %reccoef saved /; %rect saved /; call lts_rec2( ); is a subroutine version of this command /; /; example of use; /; /; b34sexec matrix; /; call loaddata; /; call load(lts :staging); /; call load(lts_rec :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); /; oldcoef=1; /; n_recur=6; /; call lts_rec; /; b34srun; /; /; Built 23 June 2007. Inprovements 1 September 2007 /; ----------------------------------------------------------- /; holdname=%names; holdlag =%lag; holdcoef=%coef; holdse =%se; holdt =%t; %reccoef=array(n_recur,nocols(newx):); %rect =array(n_recur,nocols(newx):); do i=1,n_recur; /; logic . Rerun to get %coef etc and proceed if(i.eq.1)call olsq(newy,newx,:noint :savex :holdout ihold); yy=newy; xx=newx; n_y=norows(yy); call deleterow(yy,n_y-ihold+1,ihold); call deleterow(xx,n_y-ihold+1,ihold); if(oldcoef.eq.1)then; %coef=holdcoef; %se =holdse; %t =holdt; endif; if(iprint.ne.0)then; call print(' ':); if(oldcoef.eq.0)call print('Oldcoef is Prior LTS Coef':); if(oldcoef.eq.1)call print('Oldcoef is OLS Coef':); call print('Recursive Trimmed Squares pass ',i:); endif; call lts(yy,xx,%res,%names,%lag,%coef,%se,%t,p,newy,newx,ihold,iprint); /; logic . Rerun to get %coef etc and proceed call olsq(newy,newx,:noint :savex :holdout ihold); %reccoef(i,)=%coef; %rect(i,) =%t; enddo; return; end; == ==LTS_REC2 Recursive Least Trimmed Squares - Subroutine subroutine lts_rec2(oldcoef,n_recur,p,newy,newx, holdname,holdlag,holdcoef,holdse,holdt,ihold,iprint); /; /; Does Recursive lts. Subroutine lts must be called first /; call lts_rec; is the program version of the same command that /; requires exact names for call lts( ); /; /; Needs the following set: /; oldcoef = 0 => use Base Coefficients for table /; = 1 => use prior LTS coef for table /; n_recur = => Sets number of recursive LTS estimates /; p, newy and newx set from prior call lts. This is sorted data. /; Warning: newy,newx and ihold are changed!!!!!!!! /; The next 5 arguments set base coefficients /; holdname=%names; /; holdlag =%lag; /; holdcoef=%coef; /; holdse =%se; /; holdt =%t; /; ihold => from lts /; /; iprint => from lts /; /; example of use; /; /; b34sexec matrix; /; call loaddata; /; call load(lts :staging); /; call load(lts_rec :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); /; oldcoef=1; /; n_recur=6; /; call lts_rec2(oldcoef,n_recur,p,newy,newx,%names,%lag,%coef,%se,%t, /; ihold,iprint); /; b34srun; /; /; Built 23 June 2007 as a Program /; Converted to a Subroutine 1 Septeber 2007 and renamed /; ----------------------------------------------------------- /; do i=1,n_recur; /; /; logic . Rerun to get %coef etc and proceed /; call olsq(newy,newx :noint :savex :holdout ihold); yy=newy; xx=newx; n_y=norows(yy); call deleterow(yy,n_y-ihold+1,ihold); call deleterow(xx,n_y-ihold+1,ihold); if(oldcoef.eq.1)then; %coef=holdcoef; %se =holdse; %t =holdt; endif; call print(' ':); if(oldcoef.eq.0)call print('Oldcoef is Prior LTS Coef':); if(oldcoef.eq.1)call print('Oldcoef is OLS Coef':); call print('Recursive Trimmed Squares pass ',i:); call lts(yy,xx,%res,%names,%lag,%coef,%se,%t,p,newy,newx,ihold,iprint); enddo; return; end; == ==MARS_P Display a Mars Probit Model program mars_p; /; /; Does Probit analysis and displays coef in the MARS Model /; * Adjust coef, se, a t-vals under Probit assumptions /; * -------------------------------------------------- call probit(yvar %probitx :noint :print); /; Example: /; 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; /; /; coefnew=%coef; senew=%se; tnew=%t; /; * Save adjusted probit information into Marspline model /; * ----------------------------------------------------- call restore(:file 'marss.psv'); %coef=coefnew; %se=senew; %t=tnew; call checkpoint(:file 'marszero.psv' :var %BESTIN %FLAG %DIR %CUT %YVAR %NAMES %TYPEVAR %LAG %COEF %T %MINVAR %MAXVAR %K %NOB %RSS %SUMRE %REMAX %RESVAR %MARS_VR %SE %MODTYPE %NK); /; * Display the final MarsProbit model /; * ------------------------------------------------------------------- call marspline(yvar argument(cc) :getmodel 'marszero.psv' :dispmars); return; end; == ==MARSDIAG Row Sums of Marspline Model subroutine marsdiag(xx,c_sums,r_sums,iprint,iplot,fname); /; /; Calculates row sums and col sums from a MARSPLINE model /; :xx option required to obtain %xx /; /; r_sum indicates how many vectors are active by observation /; c_sum indicate how many times a vector is non zero over all obs /; /; assume k variables on the right /; /; xx - Set as %xx if :xx option supplied to marspline /; c_sums - k element vector indicating # of times each vector ne 0 /; r_sums - n element vector indicating # of non zero vectors in row /; iprint - => 1 print c_sums r_sums /; iprint - => 2 print c_sums r_sums as columns /; iprint - => 3 print c_sums r_sums as numbered rows /; iplot - => 1 Draw a histogram /; iplot - => 2 Draw a histogram and save in file fname /; fname - => set file name to somename.wmf or somename.hp1 /; ------------------------------------------------------------------ /; test=transpose(array(nocols(xx),norows(xx):(xx.gt.0.0))); c_sums=sumcols(test); r_sums=sumrows(test); if(kind(fname).eq.-8)fname1=c1array(8:fname); if(kind(fname).eq.-1)fname1=fname; if(iprint.eq.1)then; call print('Column sums',c_sums); call print('Row sums', r_sums); endif; if(iprint.eq.2)then; call printvascmat; call print('Column sums',c_sums ); call print('Row sums',r_sums ); call printvasv; endif; if(iprint.eq.3)then; call printvasrmat; call print('Column sums',c_sums ); call print('Row sums',r_sums ); call printvasv; endif; if(iplot.eq.1)then; call graph(r_sums :plottype hist2d :pgborder :heading 'Number of non zero vectors in each observation'); endif; if(iplot.eq.2)then; ihp=0; call ialen(fname1,jjunk); ijunk=jjunk-2; jjjunk=integers(ijunk,jjunk); testname=c1array(3:); call character(t1,'hp1'); iijunk=integers(1,3); testname(iijunk)=fname1(jjjunk); call ilower(testname); if(testname(1).eq.t1(1) .and. testname(2).eq.t1(2) .and. testname(3).eq.t1(3))ihp=1; if(ihp.eq.0) call graph(r_sums :plottype hist2d :file fname1 :hardcopyfmt wmf :pgborder :heading 'Number of non zero vectors in each observation'); if(ihp.ne.0) call graph(r_sums :plottype hist2d :file fname1 :hardcopyfmt HP_GL2 :pgborder :heading 'Number of non zero vectors in each observation'); endif; return; end; == ==MARSVPLT Plots MARSSPLINE Vectors subroutine marsvplt(xx,prefix); /; Plot the individual variables in Mars Model /; xx => Obtained if :xx option used /; prefix => Prefix for the save file /; Routine built 17 July 2008 by Houston H. Stokes /; /; Note: By the nature of MARS not all vectors are no 0 all the time /; icount=0; do i=2,nocols(xx); icount=icount+1; hh1='c_______.wmf'; prefix2=c1array(8:prefix); call ialen(prefix2,jj); hh=c1array(8:prefix); jj1=integers(jj); hh1(jj1)=hh(jj1); call inttostr(icount,hh,'(i8)'); call ijuststr(hh,left); call ialen(hh,jjj); jj2=integers(jj+1,jj+jjj); jj3=integers(jjj); hh1(jj2)=hh(jj3); vector=xx(,i); title1='Mars Vector # '; call inttostr(i,hh,'(i8)'); call ijuststr(hh,left); title1=catrow(title1,hh); call graph(vector :heading title1 :nocontact :pspaceon :pgyscaleright 'ni' :pgborder :grid :nokey :plottype hist2d :pgxscaletop ' i' :colors black red :file hh1 ); icount=icount+1; enddo; return; end; == ==MARSINFO Estimate of MARSPLINE Model GCV Info program marsinfo; /; /; Calculate Final MARSPLINE GCV Info /; /; Program is experimental!! /; /; Example: /; /; b34sexec options ginclude('b34sdata.mac') member(gas); /; b34srun; /; b34sexec matrix; /; call loaddata; /; call load(marsnfo :staging); /; call echooff; /; /; call marspline(gasout gasin{1 to 6} gasout{1 to 6} :print /; :mi 6 :nk 35); /; /; call marsinfo; /; b34srun; /; +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /; call print(' ':); call print( 'Analysis of GCV, RSS and KNOT by Variable before prune step':); _gcv =%fwdinfo(,1); _rss =%fwdinfo(,2); _knot=%fwdinfo(,3); _var =%names(%iwdinfo(,1)); _lag =%lag(%iwdinfo(,1)); call tabulate(_gcv,_rss,_knot,_var,_lag); call print(' ':); call free(_n,_nn,_i_imove,_zero,_miss,_rjunk,_gcv,_rss,_knot,_var,_lag); return; end; == ==MCOVF Function for Covariance. See built in mcov( ) function mcovf(%x,%res,%lag,%damp,ifsquare); /; /; Function implements Rats command of the same name. For further /; detail see Rats Manual. /; /; Note: This function is not fast but illustrates calculation. /; For a faster implementation see built-in command mcov( ) /; that has the same arguments. /; /; mcov calculates X'uu'X with an added damping factor /; This is T*S(w) from Hansen (1982) /; /; It is assumed that olsq has been called with :savex argument to /; generate %x and %res /; /; Example: /; /; 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)); /; /; Arguments /; /; %x => x matrix set by OLSQ /; %res => Set by OLSQ. if 0.0 passed assumes uu'=1 /; %lag => Variable Lags. Set by user /; %damp => Causes expression to be multiplied by /; [(L+1-abs(k))/(L+1)]**damp where L is mag lag /; ifsquare => =0 uses u'u. ne 0 uses u /; /; Built 8 October 2005 by Houston H. Stokes /; /; ------------------------------------------------------------ k=nocols(%x); n=norows(%x); s0=missing(); if(%lag.lt.0)then; call epprint('ERROR: LAG in call mcov must be ge 0'); go to done; endif; if(%lag.lt.0.and.ifsquare.ne.0)then; call epprint('ERROR: If LAG ne 0 then ifsquare not allowed'); go to done; endif; s0=kindas(%x,matrix(k,k:)); /; notes on sx=submatrix(x,1,3,2,5); /; forms a new matrix sx containing rows 1 to 3 cols 2 to 5 /; No residual cases +++++++++++++++++++++++++++++++++++++++++++++ if(norows(%res).eq.1)then; s0=transpose(%x)*%x; if(%lag.gt.0)then; xjunk=dfloat(%lag+1); do j=1,%lag; base=kindas(%x,1.0); if(%damp.ne.0.0)then; base=kindas(%x,((xjunk-dabs(dfloat(j)))/(xjunk))**%damp); endif; /$ xt is a col do it=j+1,n; xt =transpose(submatrix(%x,it, it, 1,k)); xt_j=transpose(submatrix(%x,it-j,it-j,1,k)); /; lags > 0 lags < 0 s0 =s0 +( ((xt*transpose(xt_j))+(xt_j*transpose(xt)))*base); enddo; enddo; endif; go to done; endif; /; End of no residual cases +++++++++++++++++++++++++++++++++++++++ if(norows(%res).ne.norows(%x))then; call epprint('ERROR: norows(%x).ne.norows(%res)'); go to done; endif; if(ifsquare.eq.0)then; do it=1,n; xt =transpose(submatrix(%x,it, it, 1,k)); s0=s0+(sfam(%res(it))*sfam(%res(it))*(xt*transpose(xt))); enddo; endif; if(ifsquare.ne.0)then; do it=1,n; xt =transpose(submatrix(%x,it, it, 1,k)); s0=s0+((sfam(%res(it)))*(xt*transpose(xt))); enddo; endif; /; Residual Cases where there are lags ++++++++++++++++++++++++++++ if(%lag.gt.0)then; xjunk=dfloat(%lag+1); do j=1,%lag; base=kindas(%x,1.0); if(%damp.ne.0.0)then; base=kindas(%x,((xjunk-dabs(dfloat(j)))/(xjunk))**%damp); endif; /; Note: xt is a col do it=j+1,n; xt =transpose(submatrix(%x,it, it, 1,k)); xt_j=transpose(submatrix(%x,it-j,it-j,1,k)); /; lags > 0 lags < 0 s0=s0 +( ((xt*transpose(xt_j))+(xt_j*transpose(xt))) *base* %res(it)*%res(it-j)); enddo; enddo; endif; done continue; return(s0); end; == ==NORM NORM(X,i) of a matrix or vector function norm(x,i); /; Norm of matrix or vector /; /; matrix /; norm(x,1) largest norm1 col sum /; norm(x,2) largest singular value of x /; norm(x,3) largest norm1 row sum /; norm(x,4) Frobenius norm sqrt(sum(diag(x'*x))) /; /; vector /; norm(v,i) sum(abs(v)**i)**(1./i) i=1,.... /; norm(v,-1) max(abs(v)) /; norm(v,-2) min(abs(v)) /; /; Built 4 April 2007 /; /; This tracks 100% matlab /; /; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ kk=klass(x); if(kind(i).ne.-4)then; call epprint('ERROR: norm(x,i) must pass i as an integer'); call stop; endif; if(kk.ne.1.and.kk.ne.2)then; call epprint('ERROR: norm( ) must pass matrix or vector'); call stop; endif; if(kk.eq.2.and.i.lt.1.or.i.gt.4)then; call epprint('ERROR: norm(a,i) for matrix a must set i=1,...,4'); call stop; endif; if(kk.eq.2)then; if(i.eq.1)ans=max(sumcols(abs(x))); if(i.eq.2)ans=max(svd(x)); if(i.eq.3)ans=max(sumrows(abs(x))); if(i.eq.4)ans=sqrt(sum(diag(transpose(x)*x))); endif; if(kk.eq.1)then; if(i.eq.0 )ans=kindas(x,.1d+32); if(i.eq.-1)ans=max(abs(x)); if(i.eq.-2)ans=min(abs(x)); if(i.ge.1)then; p=kindas(x,dfloat(i)); ans=(sum((afam(abs(x))**p)))**(kindas(x,1.)/p); endif; endif; return(ans); end; == ==NTOKIN Counts number of tokins in a string subroutine ntokin(cc,nfind,itest,ibad); /; /; Counts number of tokins in a string /; /; cc => Character*1 string /; nfind => number of tokins found /; itest => 0 do not test /; 1 test for integer (5) or real (6) /; ibad => 0 all ok /; 1 problem /; /; Written 1 August 2004 by Houston H. Stokes /; Useful in reading character data into an array /; /; ++++++++++++++++++++++++++++++++++++++++++++++++++ ibase=1; nfind=0; ibad=0; if(kind(cc).ne.-1)then; call epprint('ERROR: Argument # 1 of ntokin is not character*1') ibad=1; go to finish; endif; do j=1,norows(cc)*nocols(cc); imax=0; call ib34s11(cc,ibase,ifbase,isize,itokty,inewp,imax); /; call print('isize was ',isize); /; call print('ibase was ',ibase); /; call print('ifbase was ',ifbase); /; call print('inewp was ',inewp); if(isize.eq.0.or.inewp.eq.0)go to finish; if(itest.ne. 0 .and. (dabs(itokty).ne.5.and.dabs(itokty).ne.6))then; call epprint('ERROR: A tokin is not real*8 or integer') call print( ' Type of token found ',itokty :line); i=integers(ifbase,ifbase+isize-1); find=cc(i); call character(tt,' Token found was: '); call expand(tt,find,27,(27+isize)); call print(tt :); call print(' ' :); ibad=1; go to finish; endif; i=integers(ifbase,ifbase+isize-1); ibase=inewp; nfind=nfind+1; if(inewp.eq.-99)go to finish; enddo; finish continue; return; end; == ==NW_SE Newey West SE Correction Calculation subroutine nw_se(%names,%lag,%coef,%xpxinv,%res, damp,%x,%ols_se,%w_se,%nw_se,lagnw,nwcov,whitecov,iprint); /; /; Implements Newey-West Correction to SE of a OLS Model /; In Lag = 0 get White correction /; /; Must be called after call olsq( :savex); /; Must specify lag, damp & iprint. /; damp=1.0 is rthe usual Newey-West Test /; /; For Detail see Greene Edition 5 page 267 and the Rats User's Guide /; for discussion of damp. /; /; The program rnw_se will call nw_se and set the usual defaults. /; /; It is assumed that olsq has been called with :savex argument /; /; Arguments /; /; %names => Variable Names. Set by OLSQ /; %lag => Variable Lags. Set by OLSQ /; %coef => Coef. Set by OLSQ /; %xpxinv => inv(X'X) Set by OLSQ /; %res => Residuals. Set by OLSQ /; damp => Adjusts Newey-West Weights. See Rats. /; Usual setting 1.0 /; %x => Data matrix. Set by OLSQ /; %ols_se => OLS SE Set by OLSQ /; %w_se => White SE /; %nw_se => Newey West SE /; lagnw => Lag for Newey West /; nwcov => Newey-West Var-Covar /; whitecov => White Var-Covar /; iprint => =1 iprint, =2 iprint nwcov and whitecov also k=nocols(%x); n=norows(%x); s0=kindas(%x,matrix(k,k:)); if(lagnw.le.0)then; call epprint('ERROR: LAGNW in call NW_SE must be > 0'); go to done; endif; if(nocols(%x).ne.norows(%coef))then; call epprint('ERROR: # cols of %x NE # rows %coef'); go to done; endif; /$ s0 used in White SE. For details see Greene (ed 5 eq 10-13) do it=1,n; xt =transpose(submatrix(%x,it, it, 1,k)); s0=s0+(sfam(%res(it))*sfam(%res(it))*(xt*transpose(xt))); enddo; s0=s0*kindas(%x,(1.0/dfloat(n))); if(lagnw.gt.0)then; nwcov =kindas(%x,matrix(norows(s0),nocols(s0):)); do j=1,lagnw; base=kindas(%x,(1.0 - (dfloat(j)/(dfloat(lagnw)+1.)))**damp); /$ xt is a col do it=j+1,n; xt =transpose(submatrix(%x,it, it, 1,k)); xt_j=transpose(submatrix(%x,it-j,it-j,1,k)); nwcov =nwcov +( ((xt*transpose(xt_j))+(xt_j*transpose(xt))) *base* %res(it)*%res(it-j)); enddo; enddo; nwcov=s0+(kindas(%x,(1.0/dfloat(n)))*nwcov); nwcov =kindas(%x,(dfloat(n)))*%xpxinv*nwcov*%xpxinv; whitecov=kindas(%x,(dfloat(n)))*%xpxinv*s0*%xpxinv; %w_se=dsqrt(diag(whitecov)); %nw_se =dsqrt(diag(nwcov)); if(iprint.ne.0)then; %ols_t =vfam(afam(%coef)/afam(%ols_se)); %w_t =vfam(afam(%coef)/afam(%w_se )); %nw_t =vfam(afam(%coef)/afam(%nw_se )); call print(' ' :); call print('Newey - West SE calculated using Lag ',lagnw:); call print('Damp was set as ',damp:); call tabulate(%names,%lag,%coef,%ols_se,%ols_t, %w_se,%w_t,%nw_se,%nw_t :cname,:title 'Alternative SE and t scores'); if(iprint.eq.2)call print(nwcov,whitecov); endif; endif; done continue; return; end; == ==OLS_CON Constrained OLS Model subroutine ols_con(y,x,ols_b,con_b,small_r,big_r,ols_e,con_e); /; /; Linear Constrained OLSQ /; /; y => Left Hand side input /; x => Right Hand side n by k input /; ols_b => OLS Beta k by 1 output /; con_b => Constraned Beta k by 1 output /; small_r => left hand of linear constraint input /; big_r => Right of Linear constraint input /; constraint small_r=big_r*con_b /; ols_e => OLS Error vector output /; con_e => Constrained Error vector output /; /; Built 6 August 2006 /; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /; small_r=mfam(small_r); big_r =mfam(big_r); xpxinv = inv(transpose(mfam(x))*mfam(x)); ols_b= mfam(xpxinv*transpose(mfam(x))*to_cmatrix(y)); con_b = ols_b + xpxinv*transpose(big_r)*inv(big_r*xpxinv*transpose(big_r)) * (small_r-(big_r*to_cmatrix(ols_b))); ols_e=to_vector(to_cmatrix(y)-mfam(x)*ols_b); con_e=to_vector(to_cmatrix(y)-mfam(x)*con_b); return; end; == ==POISSON Estimation of a Poisson Model program poisson; /; /; This program will run a general poisson model /; The line dmin1(x*parm,173.67) prevents dexp from getting /; too big. Limit internally in B34S is 709 /; /; Define yfact=(-1.) * dlog(fact(y)); /; syfact=sum(yfact); /; xb, y and parm /; /; func is used in maximize command /; xb = dmin1(x*parm,173.67); func1=(((-1.)*afam(dexp(xb))) +(afam(y)*afam(xb))); func=sum(func1)+syfact; call outstring(3,3,'Function'); call outdouble(36,3,func); return; end; == ==POISS_SE Specialized SE's for Poisson Models subroutine poiss_se(x,y,yhat,coef,se1,se2,se3,t1,t2,t3,iprint); /; /; x => X matrix /; y => Y vector /; yhat => predicted Y /; exb =dexp(xb); /; coef => parameters estimated /; se1 => Greene (2008, page 915) Standard Asymptotic SE /; se2 => Greene (2008, page 915) BHHH Asymptotic SE /; se3 => Greene (2008, page 915) Robust Asymptotic SE /; t1 => t using SE1 /; t2 => t using SE2 /; t3 => t using SE3 /; iprint => NE 0 implies print SE's /; r=y-yhat; r2=vfam(afam(r)*afam(r)); se1= sqrt(dabs(diag(inv(transpose(x)*diagmat(yhat)*x)))); se2= sqrt(dabs(diag(inv(transpose(x)*diagmat(r2 )*x)))); se3= sqrt(dabs(diag(inv(transpose(x)*diagmat(yhat)*x)* (transpose(x)*diagmat(r2 )*x)* inv(transpose(x)*diagmat(yhat)*x)))); t1=vfam(afam(coef)/afam(se1)); t2=vfam(afam(coef)/afam(se2)); t3=vfam(afam(coef)/afam(se3)); if(iprint.ne.0)then; call print('Alternative Poisson SE calculations':); call print('SE1 Greene (2008, 915) Standard Asymptotic SE ':); call print('SE2 Greene (2008, 915) BHHH Asymptotic SE ':); call print('SE3 Greene (2008, 915) Robust Asymptotic SE ':); call tabulate(coef,se1,t1,se2,t2,se3,t3); endif; return; end; == ==P_L_EST Probit & Logit Estimation subroutine p_l_est(typemod,y,x,func,coef,se,t,yhat,iprint,routine); /; /; Estimates probit or logit models /; /; Version 27 October 2004 /; Version 7 November 2004 /; /; typemod => Character*8 set as 'probit' or 'PROBIT' /; 'logit' or 'LOGIT' /; y => Right hand side 0-1 variable /; x => x matrix. Constant miust be included /; func => functional value /; coef => returned coefficients /; se => returned se /; t => returned t /; yhat => returned yhat /; iprint => Character*8 variables set as 'print' or 'PRINT' /; 'noprint' or 'NOPRINT' /; routine => 0 => use cmax2, 1=> max2 /; Note: SE's of max2 may be substantially different /; from cmax2 due to way calculation is made. /; /; Notes from IMSL: MAXF2 is based on the a safeguarded /; quadratic interpolation method to find a minimum /; point of a univariate function. Both the code and /; underlying algorithm are based on the routine ZXLSF /; written by M.J.D. Powell at the University of /; Cambridge. /; /; CMAX2 uses a quasi-Newton method and an active set /; strategy to solve minimization problems subject to /; simple bounds on the variables. /; /; In some cases this will converge very slowly. /; /; P_L_EST ses utility routines probit and logit. /; /; ------------------------------------------------------------------- mask1 = afam(y.eq.1.); mask2 = afam(y.eq.0.); call olsq(y x :noint); start=%coef; b=vfam(start); ll= array(norows(start):)-.1e+10; uu= array(norows(start):)+.1e+10; b=vfam(start); if(typemod.eq.'PROBIT '.or.typemod.eq.'probit ')then; if(iprint .eq.'PRINT '.or.iprint .eq.'print ')then; if(routine.eq.0)call cmaxf2(func :name probitm :parms b :ivalue start :maxfun 20000 :maxg 20000 :maxit 100000 :lower ll :upper uu :print); if(routine.eq.1)call maxf2(func :name probitm :parms b :ivalue start :maxfun 20000 :maxg 20000 :maxit 100000 :print); endif; if(iprint .ne.'PRINT '.and.iprint .ne.'print ')then; if(routine.eq.0)call cmaxf2(func :name probitm :parms b :ivalue start :maxfun 20000 :maxg 20000 :maxit 100000 :lower ll :upper uu ); if(routine.eq.1)call maxf2(func :name probitm :parms b :ivalue start :maxfun 20000 :maxg 20000 :maxit 100000 ); endif; yhat=probnorm(mfam(x)*vfam(%coef)); endif; if(typemod.eq.'LOGIT '.or.typemod.eq.'logit ')then; if(iprint .eq.'PRINT '.or.iprint .eq.'print ')then; if(routine.eq.0)call cmaxf2(func :name logitm :parms b :ivalue start :maxfun 20000 :maxg 20000 :maxit 100000 :lower ll :upper uu :print); if(routine.eq.1)call maxf2(func :name logitm :parms b :ivalue start :maxfun 20000 :maxg 20000 :maxit 100000 :print); endif; if(iprint .ne.'PRINT '.and.iprint .ne.'print ')then; if(routine.eq.0)call cmaxf2(func :name logitm :parms b :ivalue start :maxfun 20000 :maxg 20000 :maxit 100000 :lower ll :upper uu ); if(routine.eq.1)call maxf2(func :name logitm :parms b :ivalue start :maxfun 20000 :maxg 20000 :maxit 100000 ); endif; yhat=1.0 -(1.0/(1.0+dexp(afam(mfam(x)*vfam(%coef))))); endif; func=%func; coef=%coef; se =%se; t =%t; return; end; program probitm; xb=(x*b); func= mlsum((mask1*probnorm(afam(xb)) ) + (mask2*probnorm(afam((-1.)*xb)))); call outstring(3,3,'Function'); call outdouble(36,3,func); return; end; program logitm; xb=-1.* afam(x*b); func= mlsum((mask1* (1./(1.0+dexp(xb))))+ (mask2*(1.0 -(1./(1.0+dexp(xb)))) ) ); call outstring(3,3,'Function'); call outdouble(36,3,func); return; end; == ==PANEL_LIB Panel Subroutine Library /; /; Panel subroutines /; panel_t => transpose panel /; panel_df => Diffference panel /; pfe_1way => Driver for 1 way fixed effects /; pfe_2way => Driver for 2 way fixed effects /; panel_fe => Removes Mean and does one way fixed effects. /; panel2fe => Does one way and two way fixed effects. /; Does not require transpose. More general than panel_fe subroutine panel_t(%y,%x,n,t,%ynew,%xnew); /; /; Will rotate a balanced panel dataset to vary by n /; Assumes varies by t /; reg time /; 1 1940 /; 1 1941 /; . . /; 1 1980 /; n 1940 /; n 1941 /; . . /; n 1980 /; /; %y => n*t y vector /; %x => (n*t) by k x matrix /; n => # of regions /; t => # of time periods /; %ynew => transformed %y /; %xnew => transformed %x /; /; Built 30 October 2006 /; +++++++++++++++++++++++++++++++++++++++++++++++++++++++ %ynew=%y; %xnew=%x; i_in =integers(1,(n*t),t); i_out=integers(1,n); do i=1,t; %ynew(i_out) =%y(i_in); %xnew(i_out,)=%x(i_in,); i_in=i_in+1; i_out=i_out+n; enddo; return; end; subroutine panel_df(%y,%x,n,t,%ynew,%xnew); /; /; Difference the panel by t /; Assumes panel varies by t /; reg time /; 1 1940 /; 1 1941 /; . . /; 1 1980 /; n 1940 /; n 1941 /; . . /; n 1980 /; /; %y => n*t y vector /; %x => (n*t) by k x matrix /; n => # of regions /; t => # of time periods /; %ynew => transformed %y /; %xnew => transformed %x /; /; /; Used if expect unit root by time. /; Differences by time. /; /; Usage: Can be with or without (Assumed by Rats) constant: /; /; 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); /; /; Built 30 October 2006 /; ++++++++++++++++++++++++++++++++++++++++++++++++++++++ %ynew=array(norows(%y):) + missing(); %xnew=array(norows(%x),nocols(%x):)+ missing(); k=nocols(%x); i_in=integers(1,t); i_out=integers(2,t); do i=1,n; ydif=afam(dif(%y(i_in))); %ynew(i_out)=ydif; do j=1,k; xdif=afam(dif(to_array(%x(i_in,j)))); %xnew(i_out,j)=xdif; enddo; i_in=i_in +t; i_out=i_out+t; enddo; %ynew=goodrow(%ynew); %xnew=goodrow(%xnew); return; end; subroutine pfe_1way(%yy,%xx,%names2,n1,t1,iopt); /; /; one way fixed effects /; /; Shows Raw OLS equation and adjusted SE and t scores /; /; %yy => Left hand side of an OLS model with nn*tt obs /; %xx => Right hand side of model /; %names2 => Names of the series on right /; nn => # of individuals /; tt => # of time periods /; iopt => 0 One way by individuals /; 1 One way by time /; 2 difference fixed effects /; /; Data assumed to be in form: /; /; Assumes panel varies by tt /; reg time /; 1 1940 /; 1 1941 /; . . /; 1 1980 /; n1 1940 /; n1 1941 /; . . /; n1 1980 /; /; Usage: /; /; call load(panel_lib :staging); /; call olsq(invest f c :print :savex); /; call echooff; /; call pfe_1way(%y,%x,%names,10,20,0); /; call pfe_1way(%y,%x,%names,10,20,1); /; call pfe_1way(%y,%x,%names,10,20,2); /; /; Built 1 Jan 2007 /; /; +++++++++++++++++++++++++++++++++++++++++++++++++++++++ /; %h_names=%names2(integers(1,nocols(%xx)-1)); call print(' ':); call print('+++++++++++++++++++++++++++++++++++++++++++++++':); call print(' ':); if(iopt.eq.0.or.iopt.eq.1)then; call panel2fe(%yy,%xx,n1,t1,%ynew,%xnew,iopt); if(iopt.eq.0) call print('One Way Fixed Effects Estimator by Individuals':); if(iopt.eq.1) call print('One Way Fixed Effects Estimator by time/regions':); call deletecol(%xnew,nocols(%xnew)); call olsq(%ynew %xnew :print :noint); call print(' ':); call print('Raw SE and t tests and Adjusted SE and t tests':); nt=dfloat(n1*t1); nn=nt/dfloat(t1); kk=dfloat(nocols(%xx)-1); if(iopt.eq.0)adj=(nt-kk)/((nn*dfloat(t1-1))-kk); if(iopt.eq.1)adj=(nt-kk)/((dfloat(t1)*dfloat(n1-1))-kk); %adj_se=dsqrt(adj*afam(%se)*afam(%se)); %adj_t=afam(%coef)/%adj_se; call tabulate(%h_names,%coef,%se,%t,%adj_se,%adj_t); call print(' ':); endif; if(iopt.eq.2)then; call print('Differencing of a n by t panel':); call panel_df(%yy,%xx,n1,t1,%ynew,%xnew); call print('FD Fixed Effects Estimator':); call deletecol(%xnew,nocols(%xnew)); call print(' ':); call print('Estimated with constant':); call olsq(%ynew %xnew :print); call print(' ':); call print('Estimated without constant':); call olsq(%ynew %xnew :print :noint); endif; return; end; subroutine pfe_2way(%yy,%xx,%names2,n1,t1); /; /; Two way fixed effects /; /; Shows Raw OLS equation and adjusted SE and t scores /; /; %yy => Left hand side of an OLS model with nn*tt obs /; %xx => Right hand side of model /; %names2 => Names of the series on right /; n1 => # of individuals /; t1 => # of time periods /; /; Data assumed to be in form: /; /; Assumes panel varies by t1 /; reg time /; 1 1940 /; 1 1941 /; . . /; 1 1980 /; n1 1940 /; n1 1941 /; . . /; n1 1980 /; /; Usage: /; /; call load(panel_lib :staging); /; call olsq(invest f c :print :savex); /; call echooff; /; call pfe_2way(%y,%x,%names,10,20); /; Built 1 Jan 2007 /; /; +++++++++++++++++++++++++++++++++++++++++++++++++++++++ /; %h_names=%names2(integers(1,nocols(%xx)-1)); call print(' ':); call print('+++++++++++++++++++++++++++++++++++++++++++++++':); call print(' ':); call panel2fe(%yy,%xx,n1,t1,%ynew,%xnew,2); call print('Two Way Fixed Effects Estimator panel2fe=2':); call deletecol(%xnew,nocols(%xnew)); call olsq(%ynew %xnew :print :noint); call print(' ':); call print('Raw SE and t tests and Adjusted SE and t tests':); nt=dfloat(n1*t1); nn=nt/dfloat(t1); kk=dfloat(nocols(%xx)-1); adj=(nt-kk)/(((nn-1.)*dfloat(t1-1))-kk); %adj_se=dsqrt(adj*afam(%se)*afam(%se)); %adj_t=afam(%coef)/%adj_se; call tabulate(%h_names,%coef,%se,%t,%adj_se,%adj_t); call print(' ':); return; end; subroutine panel_fe(%y,%x,n,t,%ynew,%xnew); /; /; Removes the mean from a panel. /; Resulting matrix can be used to getfixed effect model /; Assumes panel varies by t /; reg time /; 1 1940 /; 1 1941 /; . . /; 1 1980 /; n 1940 /; n 1941 /; . . /; n 1980 /; /; %y => n*t y vector /; %x => (n*t) by k x matrix /; n => # of regions /; t => # of time periods /; %ynew => transformed %y /; %xnew => transformed %x /; /; Built 30 October 2006 /; ++++++++++++++++++++++++++++++++++++++++++++++++++++++ %ynew=array(norows(%y):) ; %xnew=array(norows(%x),nocols(%x):); k=nocols(%x); i_in=integers(1,t); do i=1,n; ymean=afam(mean(%y(i_in))); %ynew(i_in)=afam(%y(i_in))-ymean; do j=1,k; xmean=afam(mean((%x(i_in,j)))); %xnew(i_in,j)=afam(%x(i_in,j))-xmean; enddo; i_in=i_in +t; enddo; return; end; subroutine panel2fe(%y,%x,n,t,%ynew,%xnew,iopt); /; /; Removes the means from a rectangular panel. /; Resulting matrix can be used to get fixed effect model /; panel_fe does individual effects a bit faster. /; panel2fe does two way fixed effects on balanced design. One /; way effects do not require transpose /; /; Assumes panel varies by t /; reg time /; 1 1940 /; 1 1941 /; . . /; 1 1980 /; n 1940 /; n 1941 /; . . /; n 1980 /; /; %y => n*t y vector /; %x => (n*t) by k x matrix /; n => # of regions /; t => # of time periods /; %ynew => transformed %y /; %xnew => transformed %x /; iopt => = 0 individual effect /; iopt => = 1 time effect /; iopt => = 2 both effects /; /; Built 5 November 2006 /; ++++++++++++++++++++++++++++++++++++++++++++++++++++++ %ynew=array(norows(%y):) ; %xnew=array(norows(%x),nocols(%x):); k=nocols(%x); if(iopt.eq.0.or.iopt.eq.2)then; i_in=integers(1,t); /; Individual effect do i=1,n; ymean=afam(mean(%y(i_in))); %ynew(i_in)=afam(%y(i_in))-ymean; /; if(iopt.eq.2) %ynew(i_in)=%ynew(i_in)+y_b_mean; do j=1,k; xmean=afam(mean((%x(i_in,j)))); %xnew(i_in,j) =to_array( %x(i_in,j))-xmean; /; if(iopt.eq.2)%xnew(i_in,j)=to_array(%xnew(i_in,j))+x_b_mean(j); enddo; i_in=i_in + t; enddo; /; do time part starting from what has been built if(iopt.eq.2)then; i_in=integers(1,n*t,t); do i=1,t; ymean =afam(mean(%y(i_in))); %ynew(i_in)=%ynew(i_in)-ymean; do j=1,k; xmean =afam(mean((%x(i_in,j)))); %xnew(i_in,j)=%xnew(i_in,j)-xmean; enddo; i_in=i_in + 1; enddo; x_b_mean=array(k:); y_b_mean=mean(%y); do j=1,k; x_b_mean(j)=mean(%x(,j)); enddo; %ynew=%ynew+y_b_mean; do j=1,k; %xnew(,j)=%xnew(,j)+x_b_mean(j); enddo; endif; endif; /; Time effect if(iopt.eq.1)then; i_in=integers(1,n*t,t); do i=1,t; ymean =afam(mean(%y(i_in))); %ynew(i_in)=to_array(%y(i_in))-ymean; do j=1,k; xmean =afam(mean((%x(i_in,j)))); %xnew(i_in,j)=to_array(%x(i_in,j))-xmean; enddo; i_in=i_in + 1; enddo; endif; return; end; == ==QR_SMALL Break Appart QR Factorization subroutine qr_small(x,q,r,q1,q2,rsmall); /; /; Break appart r and q in case where the command /; /; r=qr(x,q:); /; /; was used. /; /; x => original m by n matrix. # nows must be > # cols /; q => full q from command q=qr(x,q:); /; q is m by m /; r => full r. r is m by n /; q1 => economy q /; q2 => Rest of q /; rsmall => n by n /; /; Built 18 August 2006 /; if(norows(x).le.nocols(x))then; call epprint('ERROR: # of rows of x must be > # cols of x'); go to finish; endif; if(norows(q).ne.norows(x).or.norows(q).ne.nocols(q))then; call epprint('ERROR: # of rows, cols of q must be = # rows x'); go to finish; endif; if(norows(r).ne.norows(x).or.nocols(x).ne.nocols(r))then; call epprint('ERROR: # of rows of r ne # rows of x'); go to finish; endif; 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)); finish continue; return; end; == ==RIDGE Ridge Regression subroutine ridge(y,x,lamda,ols_c,_name,_lag,ridge_c,iscale); /; /; Ridge Regression. See Hastie-Tibshirani-Friedman (2001) page 60 /; /; y => left hand side /; x => right hand side - constant at end /; lamda => Lamda for Ridge Regression /; OLS_c => OLS Coef /; _name => Usually %name /; lag => Usually %lag /; ridge_c => Ridge Coef /; iscale => =0 do not center X matrix /; =1 center X matrix /; /; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /; Routine built 24 March 2006 /; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ newx=mfam(x); if(iscale.ne.0)then; call deletecol(newx); b0=mean(y); do i=1,nocols(newx); newx(,i)=newx(,i)-mean(newx(,i)); enddo; endif; d=diagmat(vector(nocols(newx):))+lamda; b_ridge=inv((transpose(newx)*newx) + d)*transpose(newx)*vfam(y); if(iscale.ne.0)then; ridge_c=vector(:b_ridge,b0); call print(' ':); call print('Ridge Regression: X matrix has been centered ':); endif; if(iscale.eq.0)then; ridge_c=vector(:b_ridge); call print(' ':); call print('Ridge Regression: X matrix has not been centered ':); endif; call tabulate(_name,_lag,OLS_c,ridge_c); return; end; == ==RNW_SE Run Newey West SE Correction for OLS Model /; /; Duplicate copy of nw_SE in rnw_SE file /; program rnw_se; /; /; Runs Newey - West (1987) test if prior call olsq set :savex /; /; User can optionally set iprint, lag or damp. If these are not /; set, defaults are used. /; /; damp = 1.0 => is usual setting /; lag = 4 => is default setting /; iprint = 1 => prints OLS, White and NW SE and t;s /; = 2 => Print White and NW Var-Covariance Matrix /; /; Example: /; /; call load(nw_se :staging); /; call olsq(gasout gasin :print :savex); /; call rnw_se; /; if(abs(kind(iprint)).eq.99) iprint=1; if(abs(kind(damp)) .eq.99) damp =1.0; if(abs(kind(lagnw)) .eq.99)lagnw =4; call nw_se(%names,%lag,%coef,%xpxinv,%res, damp,%x,%se,%whitese,%nwse,lagnw,nwcov,whitecov,iprint); return; end; subroutine nw_se(%names,%lag,%coef,%xpxinv,%res, damp,%x,%ols_se,%w_se,%nw_se,lagnw,nwcov,whitecov,iprint); /; /; Implements Newey-West Correction to SE of a OLS Model /; In Lag = 0 get White correction /; /; Must be called after call olsq( :savex); /; Must specify lag, damp & iprint. /; damp=1.0 is rthe usual Newey-West Test /; /; For Detail see Greene Edition 5 page 267 and the Rats User's Guide /; for discussion of damp. /; /; The program rnw_se will call nw_se and set the usual defaults. /; /; It is assumed that olsq has been called with :savex argument /; /; Arguments /; /; %names => Variable Names. Set by OLSQ /; %lag => Variable Lags. Set by OLSQ /; %coef => Coef. Set by OLSQ /; %xpxinv => inv(X'X) Set by OLSQ /; %res => Residuals. Set by OLSQ /; damp => Adjusts Newey-West Weights. See Rats. /; Usual setting 1.0 /; %x => Data matrix. Set by OLSQ /; %ols_se => OLS SE Set by OLSQ /; %w_se => White SE /; %nw_se => Newey West SE /; lagnw => Lag for Newey West /; nwcov => Newey-West Var-Covar /; whitecov => White Var-Covar /; iprint => =1 iprint, =2 iprint nwcov and whitecov also k=nocols(%x); n=norows(%x); s0=kindas(%x,matrix(k,k:)); if(lagnw.le.0)then; call epprint('ERROR: LAGNW in call NW_SE must be > 0'); go to done; endif; if(nocols(%x).ne.norows(%coef))then; call epprint('ERROR: # cols of %x NE # rows %coef'); go to done; endif; /$ s0 used in White SE. For details see Greene (ed 5 eq 10-13) do it=1,n; xt =transpose(submatrix(%x,it, it, 1,k)); s0=s0+(sfam(%res(it))*sfam(%res(it))*(xt*transpose(xt))); enddo; s0=s0*kindas(%x,(1.0/dfloat(n))); if(lagnw.gt.0)then; nwcov =kindas(%x,matrix(norows(s0),nocols(s0):)); do j=1,lagnw; base=kindas(%x,(1.0 - (dfloat(j)/(dfloat(lagnw)+1.)))**damp); /$ xt is a col do it=j+1,n; xt =transpose(submatrix(%x,it, it, 1,k)); xt_j=transpose(submatrix(%x,it-j,it-j,1,k)); nwcov =nwcov +( ((xt*transpose(xt_j))+(xt_j*transpose(xt))) *base* %res(it)*%res(it-j)); enddo; enddo; nwcov=s0+(kindas(%x,(1.0/dfloat(n)))*nwcov); nwcov =kindas(%x,(dfloat(n)))*%xpxinv*nwcov*%xpxinv; whitecov=kindas(%x,(dfloat(n)))*%xpxinv*s0*%xpxinv; %w_se=dsqrt(diag(whitecov)); %nw_se =dsqrt(diag(nwcov)); if(iprint.ne.0)then; %ols_t =vfam(afam(%coef)/afam(%ols_se)); %w_t =vfam(afam(%coef)/afam(%w_se )); %nw_t =vfam(afam(%coef)/afam(%nw_se )); call print(' ' :); call print('Newey - West SE calculated using Lag ',lagnw:); call print('Damp was set as ',damp:); call tabulate(%names,%lag,%coef,%ols_se,%ols_t, %w_se,%w_t,%nw_se,%nw_t :cname,:title 'Alternative SE and t scores'); if(iprint.eq.2)call print(nwcov,whitecov); endif; endif; done continue; return; end; == ==SVD_OLS OLS Using SVD Approach subroutine svd_ols(y,x,u,v,s,beta_svd,se_svd,pc_coef,sigmasq,rss,ibad); /; /; Uses SVD to get OLS coef. Works for real*8 & real*16. /; For a discussion of methods, see Stokes (1997) pages 262-263 /; /; Uses Linpack code. See SVD2_OLS for LAPACK code. /; /; y => left hand side /; x => right hand side /; u => left hand side of SVD /; v => transpose of right hand side svd. /; x => u * diagmat(s)*V /; s => singular values /; beta_svd => SVD Beta Coef /; se_svd => SE of SVD beta coef /; pc_coef => PC Coef /; sigmasvd => sigmasq from SVD /; rss => Residual sum of squares /; /; Built 1 August 2004 /; s=svd(x,ibad,21,u,v :linpack); /; Uses fact that inv(transpose(x)*x) = v*(s**-2)*transpose(v) if(ibad.eq.0)then; pc_coef =transpose(u)*mfam(y); beta_svd=v*inv(diagmat(s))*pc_coef; rss=sumsq(y-x*beta_svd); sigmasq=rss/kindas(x,dfloat(norows(x)-nocols(x))); se_svd=vector(norows(v):); h=inv(diagmat(s)); invxpx= v*(h*h)*transpose(v); do j=1,norows(v); se_svd(j)=dsqrt(invxpx(j,j)*sigmasq); enddo; endif; return; end; == ==SVD2_OLS OLS Using SVD Approach and LAPACK subroutine svd2_ols(y,x,u,v,s,beta_svd,se_svd,pc_coef,sigmasq,rss,ibad); /; /; Uses SVD to get OLS coef. Works only for real*8 since it uses the /; LAPACK DGESVD routine which is more accurate than the LINPACK SVD /; routine. /; For a discussion of methods, see Stokes (1997) pages 262-263 /; /; y => left hand side /; x => right hand side /; u => left hand side of SVD /; v => transpose of right hand side svd. /; x => u * diagmat(s)*V /; s => singular values /; beta_svd => SVD Beta Coef /; se_svd => SE of SVD beta coef /; pc_coef => PC Coef /; sigmasvd => sigmasq from SVD /; rss => Residual sum of squares /; /; Built 1 September 2004 /; s=svd(x,ibad,21,u,v :lapack); /; Uses fact that inv(transpose(x)*x) = U*(s**-2)*transpose(v) if(ibad.eq.0)then; pc_coef =transpose(u)*mfam(y); beta_svd=v*inv(diagmat(s))*pc_coef; rss=sumsq(y-x*beta_svd); sigmasq=rss/dfloat(norows(x)-nocols(x)); se_svd=vector(norows(v):); h=inv(diagmat(s)); invxpx= v*(h*h)*transpose(v); do j=1,norows(v); se_svd(j)=dsqrt(invxpx(j,j)*sigmasq); enddo; endif; return; end; == ==TESTACC Test Accuracy using Residual Orthognality subroutine testacc(x,res,itest); /; /; Tests accuracy of OLS Calculation Via Correlation of residual and X /; /; The following routines must be loaded /; call load(cov :staging); /; call load(cor :staging); /; call load(norm :staging); /; call load(cond :staging); /; /; x => n by k x matrix with constant in col k. Usually %x /; res => residual from Model. Usually %res /; itest => =0 give norms and conditions as well as correlation test /; itest => =1 give only correlation test /; itest => =2 Give only residual correlations /; This routine 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 /; /; Routine built 1 June 2007. Requires :savex on call olsq( ) /; if(itest.eq.0)then; xpx=transpose(x)*x; call print(' ':); if(kind(xpx).ne.88)then; 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):); endif; if(kind(xpx).eq.88)then; call print('Norm 1 xpx ',norm(xpx,1)); 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 3 xpx ',cond(xpx,3)); call print('Cond 4 xpx ',cond(xpx,4)); endif; if(kind(xpx).eq.8)then; test=inv(xpx,testcond :gmat); testcond=1.0/testcond; call print('Condition from LAPACK ',testcond:); endif; if(kind(xpx).eq.16.or.kind(xpx).eq.4.or.kind(xpx).eq.8)then; test=inv(xpx,testcond ); testcond=kindas(testcond,1.0)/testcond; call print('Condition from LINPACK ',testcond:); if(kind(xpx).eq.8.or.kind(xpx).eq.16)then; e=seig(xpx); ee=max(e); ee2=min(e); eee=ee/ee2; call print('Eigen max / min ',eee:); endif; endif; if(kind(x).ne.88)then; 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):); endif; if(kind(x).eq.88)then; call print('Norm 1 x ',norm(x,1)); call print('Norm 3 x ',norm(x,3)); call print('Norm 4 x ',norm(x,4)); endif; endif; i=nocols(x); test=x; test(1,i)=res; call print(' ':); if(itest.eq.0.or.itest.eq.1)then; call print('Correlation between residual in last Col and X Matrix':); call print(cor(test)); endif; if(itest.eq.2)then; ans=cor(test); j=nocols(ans); acctest=ans(norows(ans,)); if(kind(ans).eq.8.or.kind(ans).eq.16)then; acctest(j)=kindas(test,missing()); acctest=goodrow(acctest); endif; call print('Correlation between x(1),...,x(k-1) and residual':); call print(acctest); endif; return; end; == ==TLOGIT Tests Logit/Probit Model subroutine tlogit(yval,yhat,upper,lower,desc,ntruer,ntruep nfalser,nfalsep,nunclear,ptruer,pfalser,iprint); /; -- Predicted -- /; yhat(=upper) yhat(unk) /; yval(Neg) A B U1 /; yval(Pos) C D U2 /; yval => Actual y value (0-1) /; yhat => Predicted y value in range 0.0 to 1.0 /; upper => User selection for true /; lower => User selection for false /; desc => Table description /; nfalser => A /; nFPos => B /; nFFals => C /; ntruer => D /; nTrues => C+D+U2 Number of yval=1 /; nFalses => A+B+U1 Number of yval=0 /; ntruep => B+D Number of yhat>=upper /; nfalsep => A+C Number of yhat U1+U2 Number of yhat unclassified /; nNFals => U1 yhat(unk) when yval=0 /; nNPos => U2 yhat(unk) when yval=1 /; pAC => (A+D)/(A+B+C+D) Accuracy /; ptruer => D/(C+D) True Positive Ratio /; pfalser => A/(A+B) True Negative Ratio /; pFPos => B/(A+B) False Positive Ratio /; pFFals => C/(C+D) False Negative Ratio /; pPrecise=> D/(B+D) Precision /; gMean1 => SQRT(ptruer*pPrecise) g-mean(1) /; gMean2 => SQRT(ptruer*pfalser) g-mean(2) /; Routine built 26 October 2004 /; Routine modified 05 May 2008 (BL) /; ********************************************************* ptruer =0.0; pfalser =0.0; pNPos =0.0; pNFals =0.0; pAC =0.0; pFPos =0.0; pFFals =0.0; pPrecise=0.0; nTrues = sum(yval .eq. 1.0); nFalses = sum(yval .eq. 0.0); ntruep = sum( yhat .gt. upper); nfalsep = sum( yhat .lt. lower); nNPos = sum((( yhat .lt. upper) .and. (yhat .gt.lower)) .and. (yval .eq. 1.0)); nNFals = sum((( yhat .lt. upper) .and. (yhat .gt.lower)) .and. (yval .eq. 0.0)); ntruer = sum(( yhat .gt. upper) .and. (yval .eq. 1.0)); nfalser = sum(( yhat .lt. lower) .and. (yval .eq. 0.0)); zchk=sum(yval .eq. 1.0) - nNPos; if(zchk.gt.0.0)then; ptruer = ntruer / (sum(yval .eq. 1.0) - nNPos); endif; zchk=sum(yval .eq. 0.0) - nNFals; if(zchk.gt.0.0)then; pfalser = nfalser / (sum(yval .eq. 0.0) - nNFals); endif; nunclear= idint(ntrues + nfalses - ntruep - nfalsep); if(ntruep.gt.0.0)then; pNPos = nNPos/ntruep; endif; if(nFalsep.gt.0.0)then; pNFals = nNFals/nfalsep; endif; zchk=dfloat(norows(yval)-nunclear); if(zchk.gt.0.0)then; pAC= (ntruer + nfalser)/dfloat(norows(yval)-nunclear); endif; nFPos = nTrues - ntruer - nNPos; if(nTrues.gt.0.0)then; pFFals = nFpos/nTrues; endif; nFFals = nFalses - nfalser - nNFals; zchk=nfalser + nFFals; if(zchk.gt.0.0)then; pFPos = nFFals/(nfalser + nFFals); endif; if(ntruep.gt.0.0)then; pPrecise= ntruer/ntruep; endif; gMean1 = dsqrt(ptruer * pPrecise); gMean2 = dsqrt(ptruer * pfalser); if(iprint.ne.0)then; call fprint(:clear); call fprint(:col 1 :string ' ' :print :clear); call fprint(:col 1 :string '-----------------------------------------------------' :print :clear); call fprint(:col 1 :string ' Confusion Matrix ' :print :clear); call ialen(desc,itmp); ipos=idint(dfloat(54-itmp)*.5); if(ipos .le.0) ipos=1; call fprint(:col ipos :string desc :print :clear); call fprint(:col 1 :string '-----------------------------------------------------' :print :clear); call fprint(:col 28 :string 'Predicted' :print); call fprint(:clear); call fprint(:col 16 :string 'Negative' :col 28 :string 'Positive' :col 40 :string 'Unclassified' :print); call fprint(:clear); call fprint(:col 1 :string 'Negative' :col 14 :display idint(nfalser ) '(I10)' :col 26 :display idint(nFFals) '(I10)' :col 38 :display idint(nNFals) '(I10)' :print); call fprint(:clear); call fprint(:col 1 :string 'Positive' :col 14 :display idint(nFPos) '(I10)' :col 26 :display idint(ntruer) '(I10)' :col 38 :display idint(nNPos) '(I10)' :print); call fprint(:clear); call fprint(:col 1 :string ' ' :print); call fprint(:clear); n=norows(yval); neff=n-nunclear; /$ call print('Total Number Cases : ',n: ); /$ call print('Total Classified Cases : ',neff:); /$ call print('Unclassified Positive Cases : ', idint(nNPos):); /$ call print('Unclassified Negative Cases : ', idint(nNFals):); /$ call print('Cut-off for True (>) : ',upper:); /$ call print('Cut-off for False (<) : ',lower:); /$ call print('Accuracy Rate : ',pAC:); /$ call print('True Positive Rate : ',ptruer:); /$ call print('False Positive Rate : ',pFPos:); /$ call print('True Negative Rate : ',pfalser:); /$ call print('False Negative Rate : ',pFFals:); /$ call print('Precision Rate : ',pPrecise:); /$ call print('g-mean1 : ',gMean1:); /$ call print('g-mean2 : ',gMean2:); call fprint(:col 1 :string 'Total Number Cases' :col 29 :string ':' :col 31 :display n '(i10)' :print :clear); call fprint(:col 1 :string 'Total Classified Cases' :col 29 :string ':' :col 31 :display neff '(i10)' :print :clear); call fprint(:col 1 :string 'Unclassified Positive Cases' :col 29 :string ':' :col 31 :display idint(nNPos) '(i10)' :print :clear); call fprint(:col 1 :string 'Unclassified Negative Cases' :col 29 :string ':' :col 31 :display idint(nNFals) '(i10)' :print :clear); call fprint(:col 1 :string 'Cut-off for True (>)' :col 29 :string ':' :col 31 :display upper '(f10.3)' :print :clear); call fprint(:col 1 :string 'Cut-off for False (<)' :col 29 :string ':' :col 31 :display lower '(f10.3)' :print :clear); call fprint(:col 1 :string 'Accuracy Rate' :col 29 :string ':' :col 31 :display pAC '(f10.3)' :print :clear); call fprint(:col 1 :string 'True Positive Rate' :col 29 :string ':' :col 31 :display ptruer '(f10.3)' :print :clear); call fprint(:col 1 :string 'False Positive Rate' :col 29 :string ':' :col 31 :display pFPos '(f10.3)' :print :clear); call fprint(:col 1 :string 'True Negative Rate' :col 29 :string ':' :col 31 :display pfalser '(f10.3)' :print :clear); call fprint(:col 1 :string 'False Negative Rate' :col 29 :string ':' :col 31 :display pFFals '(f10.3)' :print :clear); call fprint(:col 1 :string 'Precision Rate' :col 29 :string ':' :col 31 :display pPrecise '(f10.3)' :print :clear); call fprint(:col 1 :string 'g-mean1' :col 29 :string ':' :col 31 :display gMean1 '(f10.3)' :print :clear); call fprint(:col 1 :string 'g-mean2' :col 29 :string ':' :col 31 :display gMean2 '(f10.3)' :print :clear); endif; return; end; subroutine TLGTR(yval,yhat,ncut,ptruer,pfalser,pFPos,pFFals, accuracy,pprecise,gmean1,gmean2,rprob1,rprob2,up,lo,desc1, iCut,iMethod,iprint); /; /; This routine computes the ratio statistics of a confusion /; matrix based on a cut-off probability range of 0-1. /; /; INPUTS /; yval => Actual y value (0-1) /; yhat => Predicted y value in range 0.0 to 1.0 /; ncut => Number of cuts in range 0-1 /; icut => 0=use up and lo prob cut-offs provided /; 1=use rprob1 as prob cut-off /; 2=use rprob2 as prob cutoff /; iMethod => 1=MARS, 2=LOGIT, 3=PROBIT, 4=OLS, 5=GAM /; 6=ACE, 7=PPREG /; iPrint => 0=noprint/graph, 1=print, 2=print/graph, /; 3=2 and show graph /; OUTPUTS /; ptruer(ncut) => D/(C+D) True Positive Ratio /; pfalser(ncut) => A/(A+B) True Negative Ratio /; pFPos(ncut) => B/(A+B) False Positive Ratio /; pFFals(ncut) => C/(C+D) False Negative Ratio /; Accuracy(ncut) => (A+D)/(A+B+C+D) Accuracy /; pPrecise(ncut) => D/(B+D) Precision /; gMean1(ncut) => SQRT(ptruer*pPrecise) g-mean(1) /; gMean2(ncut) => SQRT(ptruer*pfalser) g-mean(2) /; rprob1 => Cutoff (average(rcut) where max(gmean1)) /; rprob2 => Cutoff (average(rcut) where max(gmean2)) /; Routine built 06 May 2008 /; Routine modified 25 June 2008 /; ********************************************************* if(iMethod.eq.1)then; call character(gFile,'CUTMARS.wmf'); call character(desc, 'MARS True/False Probability CutOff Range'); call character(descg, 'MARS Probability Cut-Point for True/False Classification'); endif; if(iMethod.eq.2)then; call character(gFile,'CUTPROB.wmf'); call character(desc, 'LOGIT True/False Probability CutOff Range'); call character(descg, 'LOGIT Probability Cut-Point for True/False Classification'); endif; if(iMethod.eq.3)then; call character(gFile,'CUTPROB.wmf'); call character(desc, 'PROBIT True/False Probability CutOff Range'); call character(descg, 'PROBIT Probability Cut-Point for True/False Classification'); endif; if(iMethod.eq.4)then; call character(gFile,'CUTPROB.wmf'); call character(desc, 'OLS True/False Probability CutOff Range'); call character(descg, 'OLS Probability Cut-Point for True/False Classification'); endif; if(iMethod.eq.5)then; call character(gFile,'CUTGAM.wmf'); call character(desc, 'GAM True/False Probability CutOff Range'); call character(descg, 'GAM Probability Cut-Point for True/False Classification'); endif; if(iMethod.eq.6)then; call character(gFile,'CUTACE.wmf'); call character(desc, 'ACE True/False Probability CutOff Range'); call character(descg, 'ACE Probability Cut-Point for True/False Classification'); endif; if(iMethod.eq.7)then; call character(gFile,'CUTPPREG.wmf'); call character(desc, 'ACE True/False Probability CutOff Range'); call character(descg, 'PPREG Probability Cut-Point for True/False Classification'); endif; rcut=array(ncut:); ptruer=array(ncut:); pfalser=array(ncut:); pprecise=array(ncut:); accuracy=array(ncut:); pFPos=array(ncut:); pFFals=array(ncut:); gmean1=array(ncut:); gmean2=array(ncut:); call print(' '); do i=1,ncut; rcut(i)=(dfloat(i)/dfloat(ncut)); call tlogit2(yval,yhat,rcut(i),rcut(i),desc,ntruer1,ntruep1, nfalser1,nfalsep1,nunclear1,ptruer1,pfalser1, accu1,precise1,pFPos1,pFFals1,gmeana,gmeanb,0); ptruer(i)=ptruer1; pfalser(i)=pfalser1; pprecise(i)=precise1; Accuracy(i)=accu1; pFFals(i)=pFFals1; pFPos(i)=pFPos1; gmean1(i)=gmeana; gmean2(i)=gmeanb; enddo; if(iprint.ne.0)then; ilng=integers(1,88); call character(dash,'-'); dashes(ilng)=dash; call print(dashes:); call ialen(desc,itmp); ipos=idint(dfloat(88-itmp)*.5); if(ipos .le.0) ipos=1; call fprint(:col ipos :string desc :print :clear); call print(dashes:); call fprint(:clear :col 1 :display ' Cut' :col 7 :display ' True-Pos' :col 19 :display ' True-Neg' :col 30 :display ' False-Pos' :col 41 :display ' False-Neg' :col 51 :display ' Accuracy' :col 61 :display ' Precision' :col 71 :display ' GMean1' :col 80 :display ' GMean2' :print); do i=1,ncut; call fprint(:clear :col 2 :display rcut(i) '(f5.3)' :col 9 :display ptruer(i) '(f8.3)' :col 21 :display pfalser(i) '(f8.3)' :col 32 :display pFPos(i) '(f8.3)' :col 43 :display pFFals(i) '(f8.3)' :col 53 :display accuracy(i) '(f8.3)' :col 63 :display pprecise(i) '(f8.3)' :col 72 :display gmean1(i) '(f8.3)' :col 81 :display gmean2(i) '(f8.3)' :print); enddo; endif; zmax=max(gmean1); imin=0; imax=0; do i=1,ncut; if(gmean1(i) .eq. zmax .and. imin.eq.0) imin=i; if(gmean1(i) .eq. zmax) imax=i; enddo; indx=integers(imin,imax); rprob1=mean(rcut(indx)); zmax=max(gmean2); imin=0; imax=0; do i=1,ncut; if(gmean2(i) .eq. zmax .and. imin.eq.0) imin=i; if(gmean2(i) .eq. zmax) imax=i; enddo; indx=integers(imin,imax); rprob2=mean(rcut(indx)); if(icut.eq.1)then; up=rprob1; lo=rprob1; endif; if(icut.eq.2)then; up=rprob2; lo=rprob2; endif; if(iprint.ne.0)then; call tlogit(yval,yhat,up,lo,desc1,tmp1,tmp2 tmp3,tmp4,tmp5,tmp6,tmp7,1); endif; if(iprint.eq.2)then; call graph(rcut, gmean1, gmean2 :noshow :xlabel 'cut' 'c3' :ylabelleft 'g_mean' 'c9' :plottype xyplot :nocontact :nolabel :markpoint 1 1 3 14 :colors black bblue red :hardcopyfmt WMF :pspaceon :pgyscaleright 'i' :pgxscaletop 'i' :pgborder :file gfile :heading descg); endif; if(iprint.eq.3)then; call graph(rcut, gmean1, gmean2 /; :noshow :xlabel 'cut' 'c3' :ylabelleft 'g_mean' 'c9' :plottype xyplot :nocontact :nolabel :nokey :markpoint 1 1 3 14 :colors black bblue red :hardcopyfmt WMF :pspaceon :pgyscaleright 'i' :pgxscaletop 'i' :pgborder :file gfile :heading descg); endif; return; end; subroutine tlogit2(yval,yhat,upper,lower,desc,ntruer,ntruep nfalser,nfalsep,nunclear,ptruer,pfalser, pAC,pprecise,pFPos,pFFals,gmean1,gmean2,iprint); /; -- Predicted -- /; yhat(=upper) yhat(unk) /; yval(Neg) A B U1 /; yval(Pos) C D U2 /; yval => Actual y value (0-1) /; yhat => Predicted y value in range 0.0 to 1.0 /; upper => User selection for true /; lower => User selection for false /; desc => Table description /; nfalser => A /; nFPos => B /; nFFals => C /; ntruer => D /; nTrues => C+D+U2 Number of yval=1 /; nFalses => A+B+U1 Number of yval=0 /; ntruep => B+D Number of yhat>=upper /; nfalsep => A+C Number of yhat U1+U2 Number of yhat unclassified /; nNFals => U1 yhat(unk) when yval=0 /; nNPos => U2 yhat(unk) when yval=1 /; pAC => (A+D)/(A+B+C+D) Accuracy /; ptruer => D/(C+D) True Positive Ratio /; pfalser => A/(A+B) True Negative Ratio /; pFFals => C/(C+D) False Negative Ratio /; pFPos => B/(A+B) False Positive Raio /; pPrecise=> D/(B+D) Precision /; gMean1 => SQRT(ptruer*pPrecise) g-mean(1) /; gMean2 => SQRT(ptruer*pfalser) g-mean(2) /; Routine built 05 May 2008 /; Routine modified /; ********************************************************* ptruer =0.0; pfalser =0.0; pNPos =0.0; pNFals =0.0; pAC =0.0; pFPos =0.0; pFFals =0.0; pPrecise=0.0; nTrues = sum(yval .eq. 1.0); nFalses = sum(yval .eq. 0.0); ntruep = sum( yhat .gt. upper); nfalsep = sum( yhat .lt. lower); nNPos = sum((( yhat .lt. upper) .and. (yhat .gt.lower)) .and. (yval .eq. 1.0)); nNFals = sum((( yhat .lt. upper) .and. (yhat .gt.lower)) .and. (yval .eq. 0.0)); ntruer = sum(( yhat .gt. upper) .and. (yval .eq. 1.0)); nfalser = sum(( yhat .lt. lower) .and. (yval .eq. 0.0)); zchk=sum(yval .eq. 1.0) - nNPos; if(zchk.gt.0.0)then; ptruer = ntruer / (sum(yval .eq. 1.0) - nNPos); endif; zchk=sum(yval .eq. 0.0) - nNFals; if(zchk.gt.0.0)then; pfalser = nfalser / (sum(yval .eq. 0.0) - nNFals); endif; nunclear= idint(ntrues + nfalses - ntruep - nfalsep); if(ntruep.gt.0.0)then; pNPos = nNPos/ntruep; endif; if(nFalsep.gt.0.0)then; pNFals = nNFals/nfalsep; endif; zchk=dfloat(norows(yval)-nunclear); if(zchk.gt.0.0)then; pAC= (ntruer + nfalser)/dfloat(norows(yval)-nunclear); endif; nFPos = nTrues - ntruer - nNPos; if(nTrues.gt.0.0)then; pFFals = nFpos/nTrues; endif; nFFals = nFalses - nfalser - nNFals; zchk=nfalser + nFFals; if(zchk.gt.0.0)then; pFPos = nFFals/(nfalser + nFFals); endif; if(ntruep.gt.0.0)then; pPrecise= ntruer/ntruep; endif; gMean1 = dsqrt(ptruer * pPrecise); gMean2 = dsqrt(ptruer * pfalser); /$ call print('Total Number Cases : ',n: ); /$ call print('Total Classified Cases : ',neff:); /$ call print('Unclassified Positive Cases : ', idint(nNPos):); /$ call print('Unclassified Negative Cases : ', idint(nNFals):); /$ call print('Cut-off for True (>) : ',upper:); /$ call print('Cut-off for False (<) : ',lower:); /$ call print('Accuracy Rate : ',pAC:); /$ call print('True Positive Rate : ',ptruer:); /$ call print('False Positive Rate : ',pFPos:); /$ call print('True Negative Rate : ',pfalser:); /$ call print('False Negative Rate : ',pFFals:); /$ call print('Precision Rate : ',pPrecise:); /$ call print('g-mean1 : ',gMean1:); /$ call print('g-mean2 : ',gMean2:); return; end; == ==TRIM Remove obs from front of a series subroutine trim(nfinal,series); /; /; nfinal => # of obs in the final series /; series => series to be trimmed /; /; obs removed from front. /; For use with time series and non time series 1d data /; /; For time series data see align and tslineup /; /; Built November 2008 /; ntest=norows(series); if(ntest.eq.nfinal)go to done; if(ntest.lt.nfinal)then; call print('# rows in series < nfinal':); go to done; endif; holdser=series; call copytime(series,holdser); i=integers(1,(ntest-nfinal)); series(i)=missing(); series=goodrow(series); call copytime(holdser,series,(ntest-nfinal)); done continue; return; end; == ==TSMOOTH Test smooth Command with Defaults 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 /; if(ismooth.eq.0)then; method=namelist(NCE NCEPT AVETD MAVE DMAVE ES DES HOLT WINTERS); arg1='%xhatobs'; arg2='%xhat'; endif; if(ismooth.ne.0)then; method=namelist(N_N N_A N_M A_N A_A A_M DA_N DA_A DA_M M_N M_A M_M DM_N DM_A DM_M); arg1='%sobs'; arg2='%s'; endif; rss =array(norows(method):); mad =array(norows(method):); mse =array(norows(method):); mape=array(norows(method):); mpe =array(norows(method):); corr=array(norows(method):); grid=c1array(8:' '); if(igrid.ne.0)grid=c1array(8:':grid '); do j=minlag,maxlag; do i=1,norows(method); /; do i=1,1; if(iprint.eq.0)call smooth(x :method argument(method(i)) :lag j); if(iprint.ne.0)call smooth(x :method argument(method(i)) :print :lag j :forecast index(5,5)); if(iprintdt.ne.0)then; call tabulate(argument(arg1), argument(arg2), %actual, %error :title method(i)); call print('sumsq(%error)',sumsq(%error)); endif; if(igraph.ne.0)then; _c1=c1array(8:method(i)); _c2=c1array(8:'_1.wmf '); _c3=c1array(8:'_2.wmf '); call ialen(_c1,_ic1); call ialen(_c2,_ic2); jjj =integers(1,_ic2); jjj2=jjj+_ic1; file1=_c1; file2=_c1; file1(jjj2)=_c2(jjj); file2(jjj2)=_c3(jjj); kk=c1array(8:'lag '); jjj2=integers(_ic1+2,_ic1+5); _c1(jjj2)=kk(integers(1,4)); _ic1=_ic1+4; jjj2=integers(_ic1+2,_ic1+5); hh=c1array(8:); call inttostr(j,hh,'(i4)'); call ijuststr(hh,left); _c1(jjj2)=hh(integers(1,4)); call graph(argument(arg2),%actual :heading _c1 :nocontact :pgborder argument(grid) :nolabel :file file1); call graph(%error :heading _c1 :nocontact :pgborder argument(grid) :nolabel :file file2); endif; jj=j; if(ismooth.eq.0)jj=jj-1; rss(i) = %rss(jj+1); mad(i) = %mad(jj+1); mse(i) = %mse(jj+1); mape(i) = %mape(jj+1); mpe(i) = %mpe(jj+1); corr(i) = %corr(jj+1); /; call print(%xhatmat); enddo; /; call print(%rss); if(sum_tab.ne.0)then; call print(' ':); call print('Summary Measures for in sample Forecast for Lag',j:); call tabulate(method, RSS MAD MSE MAPE MPE CORR); endif; enddo; return; end; == ==TS_TESTS Residual Specification tests subroutine ts_tests(nacf,res,y,x,bg_order,ireset1,ireset2,ireset3, lmorder); /; /; nacf => # of lags for Box-Pierce and Ljung-box /; 0 => do not do test /; res => Residual vector usually %res /; y => left hand side usually %y /; x => x matrix %x /; bg_order => Order of Breusch-Godfrey (1978) test /; => 0 => do not do test /; ireset1 => # of lags for Ramset (1969). /; Must be in the range 1 - Nob-3*ireset1-3*ireset2 /; ireset2 => # of eq 2 lags for Ramset (1969) /; ireset3 => # of lags for Ramset functional form test. Set ge 2 /; lmorder => Order of Engle LM ARCH test /; /; Box-Pierce, Ljung-Box and DW - Serial Correlation /; Breusch(1978) Godfrey Test - Serial Correlation /; RESET - Functional Form /; ARCH - ARCH effects /; Breusch-Pagan Test Three variants - Hetroscedasticity /; White Test Three variants - Hetroscedasticity /; /; Routine built December 2005. Arguments subject to change /; /; Note: het_test and sc_test are more specialized and deal with /; heteroscedasticity and serial correlation. ts_tests is a more /; general routine that includes specification tests such as /; RESET (1969) /; /; Requires: call load(b_g_test); /; call load(lmtest); /; /; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /; Example: /; /; 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 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 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); /; /; b34srun; /; /; yhat=y-res; if(nacf.gt.0)then; call print(' ':); call print('Box-Pierce and Ljung-Box Statistics':); acf_res = acf(res,nacf,se,pacf,mq,probmq); box_p = dfloat(norows(res))*cusumsq(acf_res); lag=integers(1,nacf); call tabulate(lag,acf_res,se,pacf,mq,probmq,box_p); endif; /; DW i=integers(2,norows(res)); %dw=sumsq(res(i)-res(i-1))/sumsq(res); call print(' ':); call print('Durbin - Watson ',%dw:); /; Breusch(1978) Godfrey Test if(bg_order.gt.0)then; call print(' ':); do iorder=1,bg_order; call B_G_test(iorder,x,res,gbtest,gbprob,1,0); enddo; endif; /; RESET if(ireset1.ge.1.and.ireset2.ge.2)then; call print(' ':); call reset(res,test,ireset1,ireset2,prob :print); endif; if(ireset3.ge.2)then; do iorder=2,ireset3; call reset69(y,yhat,x,rtest,prob,iorder,1); enddo; endif; /; ARCH if(lmorder.gt.0)then; call print(' ':); call lmtest(res,lmorder,lag,test,prob,1); endif; /; Breusch-Pagan Test Three variants call print(' ':); holdres=res; holdnob=norows(res); holdx=x; sig=sumsq(holdres)/dfloat(holdnob); adj_res=(afam(holdres)*afam(holdres))/sig; z=dfloat(integers(1,holdnob)); call olsq(adj_res z ); test=(((variance(adj_res)*dfloat(norows(adj_res)-1)))-%rss)/2.0; df=1.0; call fprint(:clear :display 'Breusch-Pagan Het. test using t ' :col 35 :display test '(g12.4)' :col 50 :display 'Chi-Sqr' :col 57 :display df '(f6.0)' :col 65 :display 'Probability' :col 77 :display chisqprob(test,1.0) '(g12.4)' :print); call olsq(adj_res holdx :noint ); test=(((variance(adj_res)*dfloat(norows(adj_res)-1)))-%rss)/2.0; df=dfloat(nocols(holdx)-1); call fprint(:clear :display 'Breusch-Pagan Het. test using X ' :col 35 :display test '(g12.4)' :col 50 :display 'Chi-Sqr' :col 57 :display df '(f6.0)' :col 65 :display 'Probability' :col 77 :display chisqprob(test,df) '(g12.4)' :print); newx=array(holdnob,nocols(holdx)*2-1:); do i=1,nocols(holdx); newx(,i)=afam(holdx(,i)); enddo; j=1; do i=nocols(holdx)+1,nocols(newx); newx(,i)=(afam(holdx(,j))*afam(holdx(,j))); j=j+1; enddo; /; call print(holdx,newx); call olsq(adj_res newx :noint :qr); test=(((variance(adj_res)*dfloat(norows(adj_res)-1)))-%rss)/2.0; df=dfloat(nocols(newx)-1); call fprint(:clear :display 'Breusch-Pagan Het. test using X**2' :col 35 :display test '(g12.4)' :col 50 :display 'Chi-Sqr' :col 57 :display df '(f6.0)' :col 65 :display 'Probability' :col 77 :display chisqprob(test,df) '(g12.4)' :print); /; White (1980) Het Test call print(' ':); adj_res=(afam(holdres)*afam(holdres)); z=dfloat(integers(1,holdnob)); call olsq(adj_res z :qr); test=dfloat(%nob)*%rsq; df=1.0; call fprint(:clear :display 'White (1980) Het. test using t ' :col 35 :display test '(g12.4)' :col 50 :display 'Chi-Sqr' :col 57 :display df '(f6.0)' :col 65 :display 'Probability' :col 77 :display chisqprob(test,1.0) '(g12.4)' :print); call olsq(adj_res holdx :noint :qr); test=dfloat(%nob)*%rsq; df=dfloat(nocols(holdx)-1); call fprint(:clear :display 'White (1980) Het. test using X ' :col 35 :display test '(g12.4)' :col 50 :display 'Chi-Sqr' :col 57 :display df '(f6.0)' :col 65 :display 'Probability' :col 77 :display chisqprob(test,df) '(g12.4)' :print); newx=array(holdnob,nocols(holdx)*2-1:); do i=1,nocols(holdx); newx(,i)=afam(holdx(,i)); enddo; j=1; do i=nocols(holdx)+1,nocols(newx); newx(,i)=(afam(holdx(,j))*afam(holdx(,j))); j=j+1; enddo; call olsq(adj_res newx :noint :qr); test=dfloat(%nob)*%rsq; df=dfloat(nocols(newx)-1); call fprint(:clear :display 'White (1980) Het. test using X**2' :col 35 :display test '(g12.4)' :col 50 :display 'Chi-Sqr' :col 57 :display df '(f6.0)' :col 65 :display 'Probability' :col 77 :display chisqprob(test,df) '(g12.4)' :print); /; +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ return; end; subroutine alt_aic(nob,k,rss); /; /; calculates AIC and SIC Using Alternative Formulas /; /; nob => # of obs /; k => # right hand variables /; rss => e'e in the original equation /; /; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /; call print(' ':); call print('Alternative AIC and SIC Calculations':); call print('Note: Modler gets log AIC and SIC ':); aic= dlog((rss)/dfloat(nob)) + (2.*dfloat(k)/dfloat(nob)); call print( 'AIC log and non log (Greene(2003) page 160 Rats (UG page 197)' aic,dexp(aic)); aic2=(rss/dfloat(nob)) * dexp(2.* dfloat(k)/dfloat(nob)); call print(' ':); call print('AIC Ramanathan (1998) page 165 ', aic2:); aic3=(rss/dfloat(nob)) * dexp(2.* dfloat(k+1)/dfloat(nob)); call print('AIC Maddala (2001) page 485 ', aic3:); sic= dlog(rss/dfloat(nob))+ (dfloat(k)*dlog(dfloat(nob))/dfloat(nob)); call print('SIC log and non log (Greene(2003) page 160) ', sic,dexp(sic)); return; end; == ==TSALIGN Align Time Series Data program tsalign; /; /; Will align time series and kill obs where there is no data /; call names(:); blank=' '; test1='%NAMES% '; test2='%NAMESL%'; /; call print(%names%); do i=1,norows(%names%); itest1=0; if(%names%(i).eq.test1)itest1=1; if(%names%(i).eq.test2)itest1=1; if(itest1.ne.0)%names%(i)=blank; if(%names%(i).ne.blank)then; if(freq(argument(%names%(i))).eq.0.0) %names%(i)=blank; endif; enddo; /; remove blanks in %names% and add a blank between names call sort(%names%:); mask=(%names%.eq.blank); idrop=idint(sum(mask)); j=integers(1,norows(%names%)-idrop); %names%=%names%(j); b1=c1array(1:); cc=c1array(norows(%names%),9:); do j=1,norows(%names%); b8=c1array(8:%names%(j)); cc(j,)=c1array(9:b8,b1); enddo; call tslineup(argument(transpose(cc))); cdate=chardate(%julian%); year=fyear(%julian%); call align(cdate,year,argument(transpose(cc))); return; end; == ==UPMATRIX Unpack a matrix subroutine upmatrix(oldmat,names,lags,newlevel,newnames); /; /; Unpack matrix or array to vectors/array /; /; /; oldmat => Old Matrix /; names => Names to use /; lags => Lags of vectors. If missing ignored /; newlevel => New Level to copy vectors. usually set as 100 /; to make available to calling routine use /; level()-1 /; 12345678 /; newnames => Array of new names. If lags then "COL____1" /; If nolags then uses names /; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /; Built 21 June 2007 by Houston H. Stokes /; ------------------------------------------------------------ /; kind_in =kind(oldmat); kind_lag =kind(lags); case1=0; if(kind_lag.eq.-99)case1=1; if(case1.eq.0.and.sum(lags).eq.0)case1=1; if(case1.eq.1)then; do i=1,nocols(oldmat); call copy(oldmat(,i),argument(names(i))); call makelocal(argument(names(i)) :level level() newlevel); enddo; newnames=names; go to done; endif; newnames=c8array(nocols(oldmat):); do i=1,nocols(oldmat); temp='COL_____'; junk=c1array(8:); ttemp=c1array(1,8:temp); if(i.lt.10) then; call inttostr(i,junk,'(i1)'); ttemp(1,8)=junk(1); endif; if(i.gt.9).and.(i.lt.100)then; call inttostr(i,junk,'(i2)'); ttemp(1,8)=junk(2); ttemp(1,7)=junk(1); endif; if(i.gt.99).and.(i.lt.1000) then; call inttostr(i,junk,'(i3)'); ttemp(1,8)=junk(3); ttemp(1,7)=junk(2); ttemp(1,6)=junk(1); endif; call copy(oldmat(,i),argument(ttemp)); call makelocal(argument(ttemp) :level level() newlevel); newnames(i)=c8array(1:ttemp); enddo; go to done; done continue; return; end; == ==UNDIF Undifference a series 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 /; n=norows(dseries); nseries=array(n+order:); if(norows(base).eq.1)then; nseries(1)=base; endif; if(norows(base).gt.1)then; j=integers(norows(base)); nseries(j)=base(j); endif; /; /; recursive -- cannot vectorize /; do i=order+1,n+order; nseries(i)=dseries(i-order)+nseries(i-order); enddo; return; end; == ==WALD Wald Test of restrictions subroutine wald(coef,xpxinv,r,q,nob,resvar,waldtest,probwald,iprint); /; /; Performs Wald Test. For a Reference see Greene (ed 4) page 275) /; /; coef => coefficients /; xpxinv => Usually %xpxinv from olsq command. /; Can be %varcov2 if :white has been set. In this case /; set resvar=1.0 /; r => restriction matrix /; q => right hand side /; Restriction is R*coef=q /; nob => # of observations in original model /; resvar => residual variance. Usually %resvar /; wald => Wald test /; probwald => F(norows(q),nob-norows(coef)) /; iprint => Print Test /; /; +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /; /; Example from Greene /; /; Tests restriction /; /; 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; /; /; -------------------------------------------------------------------- /; xx=(mfam(r)*vfam(coef))-vfam(q); xx=matrix(norows(xx),nocols(xx):xx); waldtest =afam(transpose(xx)*inv(r*xpxinv*transpose(r))*xx)/ afam(resvar*dfloat(norows(q))); probwald=fprob(waldtest,dfloat(norows(q)),dfloat(nob-norows(coef))); if(iprint.ne.0)then; call fprint(:clear :display 'Wald F(2,*) Test on Restrictions ' :col 33 :display waldtest '(g16.8)' :col 50 :display 'Chi-Squared(2) Value' :col 70 :display waldtest*waldtest '(g12.4)' :col 82 :display 'Probability' :col 96 :display probwald '(g12.4)' :print); endif; return; end; ==