program t23 c************************************************* c Hinich tests for dependence using bicorrelations and correlations c M. J. Hinich Version 4-11-2010 c================================================= c c The data used is divided into either non-overlapped frames of lw c observations or overlapped frames of lw observations with an overlap of c !o observations set as shown below in t23.cnl. c Data - File 2 c c File t23.out (3) is always opened. The output in t23.out are statistics c for the 2nd and 3rd order test statistics for each data frame in the c sample and summary statistics for the whole sample. c filename.graph (4) - statistics in a data structure for graphing using Excel c or any similar graphing or chart program c filename-frames.out (7) - significant frame statistics c filename.sigframes (11) - signficant frames stacked in time order. c The column number is the frame length but the number of rows depends upon c how many frames are statistically significant for the chosen false alarm c probability (test size) that is set in each t23-fil.cnl control file. c c================================================= c c Parameters in control file RUNT23.CNL. Use spaces between entries! c A delimiter symbol * must be placed on each line. c Place comments on the lines of runt23.cnl file after the symbol *. c Runt23.cnl allows runs a number of control files with different file c names containing different data files and parameters. c The parameters for each run are echoed in filename.out where filename c is the name of the data file used. c=================================================================== c c The lines for the control file RUNT23.CNL c c The parameters & commands are found using key words whose first letter c may be either upper or lower case. The rest of the key word must be in c lower case. The symbol ! is a wildcard for the value to be set. c Examples: !nrun below is the wildcard for the number of runs of the program c using different cnl files called by runt23.cnl. The integer to the right c of the = sign in runt23.cnl is read by the program and sets the number of c runs. c In any !.cnl called by runt23.cnl, !r is the wildcard for the resolution c bandwidth. The real number to the right of !r= is the resolution bandwidth c that is used. c c IMPORTANT: Leave at least ONE SPACE after the paramter c c ALWAYS use 1.0 rather than .0 since the program searches for the number c to the left of the period in the parameter. c c A delimiter symbol * must be placed on each line. c Use the keywords followed by an = sign and the setting. Some settings are c character strings while others are real numbers. c c IMPORTANT: These four entries are followed by the line c # * End of Runt23.cnl parameter reads c c The program searches for the symbol # to delimite the parameter reads type c to the left of the delimiter *. c Comments on the lines of the cnl file are placed to the right of the c symbol *. These are not read by the program. c c=================================================================== c c Number of data runs=!nrun (!nrun - integer wildcard) c c There has to be !nrun cnl files in runt23.cnl c c The following options are applied to each data run: c============================ c c H exponent=!exp (!exp - real wildcard) For T4 test=!Texp c c !exp - a real number 0 < !exp < 0.5 that determines the number of lags used c in the portmentau test. I suggest !exp=0.4 c !Texp - a real number 0 < !Texp < 0.33 that determined the number of lags of c the 4th order tricovariance test. I suggest !Texp=0.3 c c============================ c c If you want to detrend the data or detrend the data transformed to growth c rates of the data by fitting a curvilinear trend using the first four c Legendre polynomials, type the following word on a line in runt23.cnl. c c detrend c c If the word detrend is not typed in runt23.cnl the data or the log c differences is NOT detrended. c c========================= c c Prob threshold=!rpv c c !rpv - probability level threshold for the t statistics for the iterations c of each OLS fit in the c Legendre polynomial detrending fit. c If a coefficient has a t value whose absolute value has a probability c Pr(abs(t)) < 1 - rpv, it is deleted in the next OLS fit. c The covariance matrix is adjusted for the subset of lagged values used. c c========================= c c These entries are followed by the line c # * End of Runt23.cnl setup reads c c============================ c c A list of the path & names of the cnl files that pass parameters to the c program. For example if !nrun=2: c \prog\t23-fil1.cnl * cnl file c \prog\t23-fil2.cnl * cnl file c c============================ c c This list of cnl file is followed by the line c # * End of Runt23.cnl control file reads c c=========================================================== c c The control files read above contains the following lines c A delimiter symbol * must be placed on each line. c Use the keywords followed by an = sign and the setting. Some settings are c character strings while others are real numbers. c c IMPORTANT: The symbol # in a line after these parameters are set. c The program searches for the symbol # to delimite the parameter reads type c to the left of the delimiter *. c Comments on the lines of the cnl file are placed to the right of the c symbol *. These are not read by the program. c c=========================================================== c c Column number=!nuc c c !nuc - column to be used when the data file has nc columns. If the data c only has one column then set !nuc=1. c c========================= c c Frame length=!lw c c !lw - frame length > 11 c If !lw=0 then the whole data set is one frame. c c========================= c c False alarm=!fa c c !fa - false alarm probability threshold for the C or H statistics. c The test statistics are uniform U(0,1) if the data (or residuals) is c pure noise (null hypothesis). Thus if fa=0.01 then 1% of the frames c have a statistic which is falsely significant - not pure noise. . c c========================= c c Overlap=!nover Resample=!nres c c !nover - the number of observations that are shifted for ovelapping the frames. c If !nover=0 then the frames will not be overlapped. c c Use a positive integer !nres > 10 for the number of resamples if you want c to bootstrap the correlation & bicorrelation test for the fa'th quantile c of the sorted statistics from the !nres resamples with replacement. c The bootstrap is applied to the residuals of an AR(np) fit of the data where c the lag parameter np = 10*p if p > 0 as set below, or np = 10. c Type 0 if you do not want to bootstrap the test statistics. c c========================= c c Ouput files for graphing=yes or no c c The syntax matters. If yes after = then two files will be opened. c The name of these files will be the file identification string on the c first line of the data file. c 1. The fn.graph file contains the C & H and other statistics for each frame. c c 2. The .frames file has statistics for the frames that have a statistically c significant C or H statistics. c The significance level is determined by the false alarm probability !fa c value set above #. c c========================= c c Significant frames=stack or no c c If stack after = the significant frames are written to fn.sigframes. c If this line is not typed in the fn.cnl file or if sifniciant frames=no then c this file will not be opened. The default is no and this line is NOT in fn.cnl. c c========================= c c If you want to trim the data at the top quantile %ptq type the character c string: c c quantile to trim back large outliers=%ptq c c All values larger than the %ptq will be trimmed down to that quantile. c c If %ptq=0.0 or if the character string is not typed in the cnl file then c the large outliers in the data are not trimmed. c c If you want to trim the data at the bottom quantile %pbq type the character c string: c c quantile to trim back small outliers=%pbq c c All values smaller than the %pbq will be trimmed up to that quantile. c c If %pbq=0.0 or if the character string is not typed in the cnl file then c the small outliers in the data are not trimmed. c c If %pbq=0.0 & %ptq=0.0 then the data is NOT trimmed. c c========================= c c If you want to transform the data to growth rates type the string: c growth rates c c If the character string 'rates' is present in the control file then the c data is transformed into growth rates - 10.[log(x(t))-log(x(t-1))]. c If the word rates is not present in the lines above # the raw data is used. c c========================= c c If you want to hard clip the data type the character string: c c hard clip c c If this character string is present in the control file then the data will c transformed into 1 if x(t)>= mean or -1 if x(t)< mean. c c========================= c c Fit an AR(!q) model to each frame c Type the character string 'AR(!q)' or 'ar(!q)' to fit an autoregressive AR(!q) c model to each data frame and then test the residuals for dependence. The c integer p typed between '(' & ')' is read by the program. c The R squared from the fit is given in one of the columns of the graph file. c The AR fit is not done if neither of the strings are present in the lines c above #. c IMPORTANT: q must be less than or equal to p (q<=p) c c========================= c c Output significant frames='yes' or 'no' c If yes then the file contains signficant frames stacked in time order. c The jth row is is the jth significant frame and the kth column is the kth c value of that frame. c c========================= c c Unit=!u (for the time of the frame start & end) c c Type the character string 'frequency' to interpret the sampling rate as a c multiple of kHz (1000 cycle/sec) or type a 1-7 character unit for the time c unit. For business and economics applications a time unit will be used such c as hour, day, or month. c NOTE: The third number on line 2 of the data file is the sampling interval if c unit is not frequency but a time unit. If the unit is frequency then the third c number is the sampling rate in kHz. c c========================= c c Data file name with path (80 characters max starting in col 1) c c========================= c c File read start=!s file read end=!e c !s - index of datum first read, !e - index of datum last read c If !e = 0, all the data is read. c c============================ c c Data file format for the first two lines (records) c c Line 1: A character string that identifies the data in thie file . c Line 2: If the unit is in frequency then line 2 using an EXAMPLE format is c rows=!nfl, columns=!nc, sampling rate=!sr, format=(f10.5,2x,a12) unit=kHz c Line 2: If the unit is a time unit such as day then line 2 is c rows=!nfl, columns=!nc, sampling interval=!sr, format=(f10.5,2x,a12) unit=day c c The data file contains nc time series in column form where the data c format is set in the second line of the data file. c Up to 50 columns can be read. c c nfl - no. of observations, sr - sampling rate in kHz if freq. is c !nfl - no. of observations (rows), sr - sampling rate in terms of kHz if the string c 'frequency' is typed in t23-fil.cnl or in terms of sr*unit if a time unit is set. c c Data format is a character*50 string including ( ). c If there is an 'a' in the format, then a character string is expected c for the date/time in a column to the left or right of the data column. c An example format is: (f7.2,2x,a7). The time mark is in a 7 character wide c field in a column two spaced to the right of a 7 column data field. c The time mark must be less than 21 characters long. c unit=time unit or kHz c c=========================================== integer run,wk real gs,rt character*700 name character*80 buf,cnl character*50 charead character*20 par character*10 su character*2 delim character dtnd logical exs c******************************************** open(9,file='T23-error.out') c Inquire if file exists inquire(file='Runt23.cnl',exist=exs) if(.not.exs) then write(6,'(/)') write(6,*)'Runt23.cnl file does not exist' write(9,*)'Runt23.cnl file does not exist' stop endif c Open control file open(1,file='Runt23.cnl',err=1,iostat=io1,status='old') c Start time call clock(gs,0,1) c------- ip=0 kt=0 ia=0 dowhile(kt<=17 .and. ia==0) kt=kt+1 c Clear buf do n=1,80 buf(n:n)=' ' enddo c Read line read(1,'(a80)',end=2,err=3,iostat=io1) buf ia=index(buf,'#') ig=index(buf,'*') if(ig==0) then write(6,*)'Place comment marker * after each line in Runt23.cnl' write(9,*)'Place comment marker * after each line in Runt23.cnl' write(6,'(a70)') buf write(9,'(a70)') buf stop endif ig=ig-1 name(ip+1:ip+ig)=buf(:ig) ip=ip+ig enddo c No. of runs par='uns' delim(1:2)='= ' run=numread(par,name,delim,ia,ier,io) c Exponent par='onent' exp=0.4 exp=rnumread(par,name,delim,ia,ier,io) if(ier>0) then write(6,'(/2x,''Error in reading the H exponent from Runt23.cnl'' *)') write(9,'(/2x,''Error in reading the H exponent from Runt23.cnl'' *)') stop endif if(exp<=0.01.or.exp>0.5) then write(6,*)'exponent = ',exp,' is out of bounds in Runt23.cnl' write(9,*)'exponent = ',exp,' is out of bounds in Runt23.cnl' stop endif c 4th order exponent par='est' te=0.3 te=rnumread(par,name,delim,ia,ier,io) if(ier>0) then write(6,'(/2x,''Error in reading T4 exponent from Runt23.cnl'')') write(9,'(/2x,''Error in reading T4 exponent from Runt23.cnl'')') stop endif if(te<=0.01.or.te>0.33) then write(6,*)'te exponent = ',te,' is out of bounds in Runt23.cnl' write(9,*)'te exponent = ',te,' is out of bounds in Runt23.cnl' stop endif c Detrend using a polynomial least squares fit dtnd='n' c Curvilinear detrending ia=index(name,'trend') if(ia>0) then dtnd='y' endif if(dtnd=='y') then c Probability threshold for the dtrend fit par='ty threshold' pv=rnumread(par,name,delim,ia,ier,6) if(ier>0) then write(6,*)'Error in read of prob threshold in cnl file' write(9,*)'Error in read of prob threshold in cnl file' stop endif c Trap pv if(pv<=0.) then write(3,*)'Prob level pv ',pv,' < 0' write(3,*)'Set pv = 0 if you do not want to detrend' write(6,*)'Prob level pv ',pv,' < 0' write(6,*)'Set pv = 0 if you do not want to dtrend' write(9,*)'Prob level pv ',pv,' < 0' write(9,*)'Set pv = 0 if you do not want to dtrend' stop endif if(pv>1) then write(3,*)'Prob level pv ',pv,' > 1' write(6,*)'Prob level pv ',pv,' > 1' write(9,*)'Prob level pv ',pv,' > 1' stop endif endif c------- c Multiple runs do nr=1,run c Read control file read(1,'(a80)',end=4,err=5,iostat=io2) cnl ig=index(cnl,'*') cnl=cnl(:ig-1) call top(cnl,exp,te,dtnd,pv,gs,run,nr,3) enddo c End of loop on runs close(9) call flush(3) stop 1 write(6,*) '**** Error on opening Runt23.cnl ',io1 stop 2 write(6,*) 'End of file on runt23.cnl reads' write(9,*) 'End of file on runt23.cnl reads' stop 3 write(6,*) 'Place the symbol # after the two top reads ',io1 write(9,*) 'Place the symbol # after the two top reads ',io1 stop 4 write(6,*) 'End of file on cnl read in runt23.cnl' write(9,*) 'End of file on cnl read in runt23.cnl' stop 5 write(6,*) '**** Error on reading name of cnl file ',nr,io2 write(6,*) cnl write(9,*) '**** Error on reading name of cnl file ',nr,io2 write(9,*) cnl stop c------- end c c=========================================== c subroutine top(cnl,exp,te,dtnd,pv,gs,run,nr,io) c*********************************************************** c Calls wind c=========================================================== parameter(mt=2,ndy=14,ny=2010) allocatable ::x,h,c,h4,a,s,k3,k4,rg,r2,ac,e,indx,t,rht,sht,v,w,y real x(:),h(:),c(:),h4(:),a(:),s(:),k3(:),k4(:),rg(:),r2(:),w(:),y *(:),e(:),v(:),rht(:),sht(:),exp,th,gs,rt real*8 ac(:) integer indx(:),q,lw,is,ie,n,np,ov,res,run,wk character*700 name character*80 buf,cnl,fn,id character*50 charead,fmt character*25 t(:) character*20 par character*7 unit character*2 cl,delim,gf character dtnd,su,mk logical exs intent(in) cnl,exp,te,dtnd,pv,gs,nr,run c******************************************** inquire(file=cnl,exist=exs) if(.not.exs) then write(6,'(/)') write(6,*) 'cnl file does not exist ',cnl write(9,'(/)') write(9,*) 'cnl file does not exist ',cnl return endif c Open control file open(10,file=cnl,err=1,iostat=io1,status='old') c------- ip=0 do k=1,17 c Clear buf do n=1,80 buf(n:n)=' ' enddo c Read line read(10,'(a80)') buf ia=index(buf,'#') if(ia>0) then exit endif ig=index(buf,'*') if(ig==0) then write(io,*)'Place comment marker * after parameters on t23.cnl' write(io,'(a70)')buf stop endif ig=ig-1 name(ip+1:ip+ig)=buf(:ig) ip=ip+ig enddo c------- c Data file name from control file par='ile name' delim(1:2)='= ' fn=charead(par,name,delim,ia,ier,3) if(ier>0) then write(6,*)'Error in reading the data file name from cnl file' write(9,*)'Error in reading the data file name from cnl file' stop endif close(2) c Inquire if file exists inquire(file=fn,exist=exs) if(.not.exs) then write(6,*)'Data file ',fn,' does not exist' write(9,*)'Data file ',fn,' does not exist' stop endif ipd=index(fn,'.') if(ipd==0) then write(6,*)'Data file ',fn,' has no extension' write(9,*)'Data file ',fn,' has no extension' stop endif c Parse for up to 4 path switches ib=index(fn,'\') if(ib>0) then ic=index(fn(ib:),'\') ib=ib+ic if(ic>0) then ic=index(fn(ib:),'\') ib=ib+ic if(ic>0) then ic=index(fn(ib:),'\') ib=ib+ic endif endif else ib=1 endif c Open data file open(2,file=fn,err=2,iostat=io2,status='old') c Open T23 output file buf=cnl ipd=index(buf,'.') buf(ipd+1:ipd+3)='out' open(3,file=buf(:ipd+3),err=3,iostat=io3) c Name of program write(3,'(7x,''Moving Frame Tests by M. J. Hinich'')') if(run>1) then write(3,'(11x,''File Name = '',a20)') buf(:ipd-1) endif c------- c Parse to open files *.graph & *.frames files par='graph files' c Graph file is fn.graph gf(1:1)='n' fmt=charead(par,name,delim,ia,ier,3) ic=index(fmt,'yes') if(ic>0) then gf(1:1)='y' endif c------- c Parse to open *.sigframes for signficant frames gf(2:2)='n' if(gf(1:1)=='y') then par='frames' fmt=charead(par,name,delim,ia,ier,3) ic=index(fmt,'stack') if(ic>0) then gf(2:2)='w' endif endif c------- c Open frames and graph files if open files for graphing=yes if(gf(1:1)=='y') then close(4) close(7) buf(ipd+1:ipd+7)='graph' open(4,file=buf(:ipd+7),err=4,iostat=io4) buf(ipd:ipd+11)='-frames.out' open(7,file=buf(:ipd+12),err=5,iostat=io5) if(gf(2:2)=='w') then buf(ipd:ipd+14)='-sigframes.out' open(11,file=buf(:ipd+15),err=6,iostat=io6) endif endif c------- c Column number par='number' ncm=numread(par,name,delim,ia,ier,io) if(ier>0) then par='Number' ncm=numread(par,name,delim,ia,ier,io) if(ier>0) then write(6,*)'Error in reading the column number from t23.cnl' write(6,*)'Error number for rnumread = ',ier write(9,*)'Error in reading the column number from t23.cnl' write(9,*)'Error number for rnumread = ',ier stop endif endif c Read frame length, false alarm probability (size) par='ength' lw=numread(par,name,delim,ia,ier,io) if(ier>0) then write(6,*)'Error in reading the frame length from t23.cnl' write(9,*)'Error in reading the frame length from t23.cnl' stop endif c False alarm threshold th=0.0 par='larm' th=rnumread(par,name,delim,ia,ier,io) if(ier>0) then write(6,*)'Error in reading the false alarm threshold' write(9,*)'Error in reading the false alarm threshold' stop endif c Trap frame length lw if(lw>0.and.lw<12) then write(3,*)'Frame length ',lw,' < 12. Increase lw ' stop elseif(lw<0) then write(3,*)'Frame lenght ',lw,' < 0' stop endif c Trap threshold if(th>1.0) then write(3,*)th,' False alarm threshold > 1' stop endif c------- c Time or frequency switch par='nit' su='t' fmt=charead(par,name,delim,ia,ier,3) if(ier>0) then write(6,'(/7x,''Error - Place unit= in the cnl file'')') write(3,'(/7x,''Error - Place unit= on T23.cnl file'')') stop endif ia=index(fmt,'quency') iunit=0 if(ia>0) then su='f' iunit=1 else su='t' c Set sampling unit unit(1:7)=fmt(1:7) endif if(su/='f' .and. su/='t') then write(6,*)'Time-frequency switch ',su,' is incorrect' write(3,*)'Time-frequency switch ',su,' is incorrect' stop endif c Frame overlap ov=0 par='lap' ov=numread(par,name,delim,ia,ier,io) c Resamples res=0 par='sample' res=numread(par,name,delim,ia,ier,io) if(ier>0) then par='samples' res=numread(par,name,delim,ia,ier,io) if(ier>0) then write(6,*)'Error ',ier,' in read of resamples' write(3,*)'Error ',ier,' in read of resamples' write(3,*)par stop endif endif c Trap res if(res<11) then res=0 endif c Data transformation switches cl='no' c------- ia=index(name,'ates') if(ia>0) then c Growth rates cl='rr' endif c------- ib=index(name,'clip') if(ib>0) then c Hard clip data cl='cl' endif c------- ic=index(name,'rim') c Trap ib & ic if(ib>0.and.ic>0) then write(6,*)'trim & hard clip swiches are set. Use only one' write(9,*)'trim & hard clip swiches are set. Use only one' stop endif if(ic>0) then c Trim tailes cl='tr' c Quantile for large outliers pcu=0.0 pcl=0.0 par='rge outliers' pcu=rnumread(par,name,delim,ip,ier,io) if(ier>0) then write(io,'(2x,''Error in read of large outliers quantile''/)') write(6,'(2x,''Error in read of large outliers quantile''/)') write(9,'(2x,''Error in read of large outliers quantile''/)') return endif c Quantile for small outliers par='all outliers' pcl=rnumread(par,name,delim,ip,ier,io) if(ier>0) then write(io,'(2x,''Error in read of small outliers quantile''/)') write(6,'(2x,''Error in read of small outliers quantile''/)') write(9,'(2x,''Error in read of small outliers quantile''/)') return endif c Trap trim levels if(pcu>=100.0.or.pcl<=0.0.and.ia==0) then cl='no' endif if(pcu>=100.0.or.pcl<=0.0.and.ia>0) then cl='rr' endif endif c------- if(ia>0.and.ic>0) then cl='tg' endif c------- c Fit an AR(q) to each frame q=0 ia=index(name,'AR') ib=index(name,'ar') if(ia>0) then par='R' delim(1:2)='()' q=numread(par,name,delim,ia,ier,io) elseif(ib>0) then par='r' delim(1:2)='()' q=numread(par,name,delim,ia,ier,io) if(ier>0) then write(6,'(/7x,''Error in read of AR lag '',i4)')q write(9,'(/7x,''Error in read of AR lag '',i4)')q write(6,*)par write(9,*)par stop endif endif c Trap q for AR(q) if(q<0) then write(6,*)'AR parameter ',q,' < 0' write(9,*)'AR parameter ',q,' < 0' stop endif if(q>lw/2.and.lw>0) then write(6,*)'AR parameter ',q,' > lw/2 ',lw/2 write(9,*)'AR parameter ',q,' > lw/2 ',lw/2 stop endif c------- delim(1:2)='= ' c File read start par='start' is=numread(par,name,delim,ia,ier,io) if(ier>0) then write(6,'(/7x,''Error in read of start datum '',i9)')is write(9,'(/7x,''Error in read of start datum '',i9)')is write(6,*)par write(9,*)par stop endif c File read end par='end' ie=numread(par,name,delim,ia,ier,io) if(ier>0) then write(6,'(/7x,''Error in read of end datum '',i9)')ie write(9,'(/7x,''Error in read of end datum '',i9)')ie write(3,*)par write(9,*)par stop endif if(is<1) then write(3,*) 'IS in CNL file is < 1 ',is stop elseif(ie/=0.and.ienc) then write(6,'(/7x,''Column '',i3,'' selected > no. columns '',i3)')nc *m,nc write(9,'(/7x,''Column '',i3,'' selected > no. columns '',i3)')nc *m,nc stop endif c------- if(ie==0) then ie=nfl endif n=ie-is+1 c Trap file length if(n>nfl) then write(6,*) 'n ',n,' > sample size ',nfl,' ',name write(3,*) 'n ',n,' > sample size ',nfl,' ',name stop endif c------- c Set lw=n if lw=0 if(lw==0) then lw=n endif np=1 c Determine no. of frames np if(ov>0) then c The frame will be moved ov at each iteration np=(n-lw-1)/ov+1 if(np>n) then np=np-1 endif elseif(ov==0) then c No overlap np=n/lw endif c Trap np if(np<1) then c Sample is the frame np=1 lw=n endif if(ov==0) then iw=lw elseif(ov>0) then c Set pointer for overlap iw=ov else write(6,*)'Overlap parameter ',ov,' is incorrect' write(3,*)'Overlap parameter ',ov,' is incorrect' stop endif c------- ng3=nint((real(lw)**exp)) nbc=(ng3-1)*ng3/2 if(nbc>=lw/2) then write(6,*)'Bicorrelations ',nbc,' > frame length/2 ',lw/2 write(9,*)'Bicorrelations ',nbc,' > frame length/2 ',lw/2 stop endif if(ng3<1) then write(6,*)'ng3 =',ng3,' < 1' write(6,*)'Increase lw or increase exp in parameter line' write(9,*)'ng3 =',ng3,' < 1' write(9,*)'Increase lw or increase exp in parameter line' stop endif c------- c Parse format line in buf to see if there is a 'a' for time marks c If so then truncate buf to data format and set mk='m' to flag times mk='a' ia=index(fmt,'a') if(ia>0) then mk='m' ia=index(fmt(ia+1:),')') ib=index(fmt(ia+1:),')') if(ib==0) then do i=ib+1,50 fmt(i:i)=' ' enddo elseif(ib>0) then ic=index(fmt(ib+1:),')') if(ic==0) then do i=ic+1,50 fmt(i:i)=' ' enddo elseif(ic>0) then ib=index(fmt(ic+1:),')') if(ib==0) then do i=ib+1,50 fmt(i:i)=' ' enddo endif endif endif endif c------- c Set pointers i1=n*nc-1 i2=np-1 i3=np*lw-1 i4=q i5=lw-1 c------- c Allocate space allocate(x(0:i1),h(0:i2),c(0:i2),h4(0:i2),a(0:i2),s(0:i2),k3(0:i2) *,k4(0:i2),rg(0:i2),r2(0:i2),sht(0:i2),rht(0:i2),ac(0:i4),e(0:i4),i *ndx(0:i3),v(0:i5),w(0:i3),y(0:i3),stat=ibad) if(ibad/=0) then write(3,*) 'Unable to allocate work space for arrays in top' write(6,*) 'Unable to allocate work space for arrays in top' write(9,*) 'Unable to allocate work space for arrays in top' stop endif if(mk=='m') then allocate(t(0:n-1),stat=ibad) elseif(mk/='m') then allocate(t(0:0),stat=ibad) endif if(ibad/=0) then write(3,*)'Unable to allocate work space for array t - ie ',ie write(6,*)'Unable to allocate work space for array t - ie ',ie write(9,*)'Unable to allocate work space for array t - ie ',ie stop endif c------- call arw(x,h,c,h4,a,s,k3,k4,rg,r2,w,e,v,ac,indx,y,sht,rht,unit,su, *ov,mk,cl,t,sr,th,id,fmt,is,ie,n,lw,np,q,dtnd,exp,gf,pcu,pcl,pv,res *,te,nr,nc,ncm,mt,ndy,ny,io) deallocate(x,h,c,h4,a,s,k3,k4,rg,r2,ac,e,indx,t,rht,sht,v,w,y) c+++++++++++++++++++++++++++++++++++++++++++ c Time and date of run & run time call head(mt,ndy,ny,gs,rt,3) c+++++++++++++++++++++++++++++++++++++++++++ return c------ 1 write(3,*) '**** Error on opening cnl file no. ',nr,io1 write(3,*)cnl write(6,*) '**** Error on opening cnl file no. ',nr,io1 write(6,*)cnl 2 write(3,*) '**** Error on opening data file ',io1,' run = ',nr write(3,'(a80)')fn return 3 write(3,*)'**** Error on opening output file ',io3 return 4 write(3,*)'**** Error on opening graph file ',io4,' run = ',nr return 5 write(3,*)'**** Error on opening frames file ',io5,' run = ',nr return 6 write(3,*)'**** Error on opening sig-frames file ',io6,' run = ', *nr return c------- end c c=========================================== c subroutine arw(x,h,c,h4,a,s,k3,k4,rg,r2,w,e,v,ac,indx,y,sht,rht,un *it,su,ov,mk,cl,t,sr,th,id,fmt,is,ie,n,lw,np,q,dtnd,exp,gf,pcu,pcl, *pv,res,te,nr,nc,ncm,mt,ndy,ny,io) c*********************************************************** c dtnd - curvalinear prewhitening, q - AR(q) for fitting each frame c=========================================================== real x(n,nc),h(np),c(np),h4(np),a(np),s(np),k3(np),k4(np),rg(np),r *2(np),w(lw*np),y(lw*np),e(q+1),v(lw),sht(np),rht(np),out(6) real*8 ac(q+1),cd(4) integer indx(np*lw),q,ov,res,lw,is,ie,n,np,nc,ncm,nbc,ntc,po,io character*80 id character*50 fmt character*25 t(*) character*24 dat character*7 unit character*2 cl,gf character dtnd,su,mk intent(in) unit,su,ov,mk,cl,sr,th,id,fmt,is,ie,n,lw,np,q,dtnd,exp, *gf,pcu,pcl,pv,res,te,nc,ncm,mt,ndy,ny,io intent(out) h,c,h4 c**************************** write(io,'(2x,57(''*''))') write(io,'(a68,'' Column '',i3)')id,ncm c Read data call datc(x,t,fmt,mk,is,ie,n,nc,2,io) c No. used nusd=np*lw c x(1,ncm) =>y do k=1,nusd y(k)=x(k,ncm) enddo c Compute statistics for the raw data or growth rates of the data call stat(nusd,y,gm,gs,g3,g4,g6,smax,smin,io) gr=smax-smin c------- if(cl=='rr'.or.cl=='tg') then c Growth rates y =>w call areturns(y,w,nusd,io) c Rates =>y do k=1,nusd y(k)=w(k) enddo c Compute statistics for the raw data or growth rates of the data call stat(nusd,y,gram,grsg,grc3,grc4,grc6,grmax,grmin,io) endif c------- po=0 if(dtnd=='y') then c Detrend call detrend(y,w,cd,in,nusd,pv,po,ram,rsd,io) call stat(nusd,w,wm,sigw,cw3,cw4,cw6,amax,amin,io) gm=wm gr=amax-amin gs=sigw g4=cw4 c Residuals w =>y do k=1,n y(k)=w(k) enddo endif c------- if(cl=='tr'.or.cl=='tg') then c Trim data => w call trim(y,w,nusd,pcu,pcl,wu,wl,io) nlc=nint(pcl/100.*float(n)) nuc=nint((1.0-pcu/100.0)*float(n)) call stat(nusd,w,am,sg,c3,c4,c6,cmax,cmin,io) endif c------- if(cl=='cl') then c Hard clip data => w call hardclp(y,w,nusd,gm,io) endif c------- if(cl=='no'.and. cl=='no') then c y=> w do k=1,nusd w(k)=y(k) enddo endif c------- c Set th threshold in tests to be used if res=0 if(res==0) then do k=1,6 out(k)=th enddo endif if(res>1) then c Resample w for test thresholds call resamp(w,v,nusd,lw,res,exp,te,gm,gs,g4,th,out,ir,io) endif c------- call wind(w,h,c,h4,a,s,k3,k4,rg,r2,w,e,v,ac,indx,unit,su,ov,mk,t,s *r,th,nusd,lw,nbc,ntc,np,q,gf,nr,exp,te,gm,gr,gs,g4,sht,rht,out,sum *fs,sumfn,sefs,sefn,nch,ncc,io) c------- c Write header for summary stats in file 3 call frwrt(fmt,unit,su,ov,cl,dtnd,sr,th,exp,is,ie,lw,np,res,out,q, *te,wl,wu,pcu,pcl,nuc,nlc,nbc,ntc,io) write(io,'(60(''_'')/)') c------- if(gf(1:1)=='y') then c Write header for frame stats file 7 write(7,'(a68,'' Column '',i3)')id,ncm call frwrt(fmt,unit,su,ov,cl,dtnd,sr,th,exp,is,ie,lw,np,res,out,q *,te,wl,wu,pcu,pcl,nuc,nlc,nbc,ntc,7) write(7,'(a1)') write(7,'('' Version '',i2,''-'',i2,''-'',i4)')mt,ndy,ny c Output time and date call fdate(dat) write(7,'('' Date & Time '',a24)')dat write(7,'(a1)') endif write(io,'(10x,''Descriptive Statistics of Data'')') call wrstat(io,gm,gs,g3,g4,g6,smax,smin) if(cl=='rr'.or.cl=='tg') then write(io,'(2x,''Statistics of the Growth Rates'')') call wrstat(io,gram,grsg,grc3,grc4,grc6,grmax,grmin) endif if(cl=='tr'.or.cl=='tg') then write(io,'(2x,''Statistics of the Trimmed Series'')') call wrstat(io,am,sg,c3,c4,c6,cmax,cmin) endif c Write header for graph file 4 if(gf(1:1)=='y') then write(4,'(a68,'' Column '',i3)')id,ncm call frwrt(fmt,unit,su,ov,cl,dtnd,sr,th,exp,is,ie,lw,np,res,out,q *,te,wl,wu,pcu,pcl,nuc,nlc,nbc,ntc,4) write(4,'(a1)') write(4,'('' Version '',i2,''-'',i2,''-'',i4)')mt,ndy,ny write(4,'('' Date & Time '',a24)') dat write(4,'(a1)') endif if(nch>9) then c Correlations of test statistics call cors(h,c,h4,s,k3,k4,rg,r2,indx,sumfs,sumfn,sefs,sefn,q,np,lw * ,nch,ncc,io) endif if(res>0) then write(3,'(/2x,''Sig-CI are standard deviations that are outside a * 99% bootstrapped confidence interval'',/'' of its bootstrapped me *an'')') write(3,'(2x,''Ran-CI are ranges that are outside a 99% bootstrap *ped confidence interval of its'',/'' bootstrapped mean'')') endif write(3,'('' The mean (Mean), standard dev (Sig), kurtosis (Kurt) *and range for each frame are'')') write(3,'(''divided by their values computed from all the data to *normalize these statistics for plotting'')') write(3,'(4x,''Mshift is the Farley-Hinich test of a mean shift'') *') write(3,'(4x,''Vshift is the Farley-Hinich test of a variance shif *t'')') write(3,'(7x,''The values are cdf probabilities under the null of *no shift''/)') c------- return end c c=========================================== c subroutine wind(x,h,c,h4,a,s,k3,k4,rg,r2,w,e,v,ac,indx,unit,su,ov, *mk,t,sr,th,n,lw,nbc,ntc,np,q,gf,nr,exp,te,gm,gr,gs,g4,sht,rht,out, *sumfs,sumfn,sefs,sefn,nch,ncc,io) c*********************************************************** c Compute arrays of stats for frames of length lw c Input: x - data array, w - frames, mk='y' - times c unit - data unit, sr - sampling rate/interval, mk - time mark c np - no. of frames, su - 't' or 'f', ar='y' - resids, c ov='y' - overlap c Output: b - bicorrelations for frame, r - correlations for frame c h - bicovariance test statistic, t4 - tricovriance test statistic c c - correlation test statistic, a - means, s- sigs, k3 - skews, c k4 - kurtoses, rg - ranges, nbc - H lags, ntc - T lags c sumfs,sefs - sum & sum of squares of absolute errors of forecasts for c signficicant H frames. sumfn,sefn - sum & sum of squares of absolute errors c of forecasts for insignficicant H frames. c=========================================================== real x(n),h(np),c(np),h4(np),a(np),s(np),k3(np),k4(np),rg(np),r2(n *p),w(lw,np),e(q+1),v(lw),sht(np),rht(np),out(6) real*8 ac(q+1),s2,s3 integer indx(np*lw),q,lw,ng3,ng4,nbc,ntc,n,np,ov,io character*25 t(*) character*7 unit character*2 gf character su,mk intent(in) x,unit,su,ov,mk,sr,th,n,lw,np,nr,q,gf,gm,gr,gs,g4,out,e *xp,te,io intent(out) h,c,h4,sumfs,sumfn,sefs,sefn,nbc,ntc,nch,ncc c**************************** c nch - no of frames with significant H statistics c ncc - no of frames with significant C statistics c ntt - no of frames with significant T4 statistics nch=0 ncc=0 nt4=0 ntt=0 if(ov==0) then iw=lw elseif(ov>0) then iw=ov else write(io,*)'Overlap parameter ',ov,' is incorrect' stop endif c------- c Pointer for last frame in w ng3=nint(float(lw)**exp) ng4=nint(float(lw)**te) ipw=(np-1)*iw iaw=((ng3-1)*ng3/2)+1 sumfs=0. sumfn=0. sefs=0. sefn=0. c Loop on np frames do kp=1,np nw=(kp-1)*iw+1 nn=nw-1 flg=1.0 ar2=0. c If q>0 compute AR residuals for frame, else normalize x => v if(q>0) then call arcs(x(nw),ac,v,e,ar2,se,ymn,lw,q,ier,io) if(ier>0) then write(io,'(2x,''Error in sub arcs frame = '',i6,'' 1st datum = *'',i7,'' run '',i2/)')kp,nn,nr flg=0.0 endif call forecast(x(nw),ac,ymn,fcy,lw,q) c One step ahead forcast error fce=abs(x(nw+1)-fcy) call stat(lw,v,am,sig,sk,c4,c6,smax,smin,io) rng=smax-smin if(sig>0.0) then c Normalize residuals do j=1,lw v(j)=(v(j)-am)/sig enddo endif if(sig<=0.0) then flg=0.0 endif endif c End of AR(q) fit if(q==0) then c No AR fit of frame call arcs(x(nw),ac,v,e,ar2,se,ymn,lw,1,ier,io) if(ier>0) then write(io,'(2x,''Error in sub arcs frame = '',i6,'' 1st datum = *'',i7,'' run '',i2/)')kp,nn,nr flg=0.0 endif call forecast(x(nw),ac,ymn,fcy,lw,1) c One step ahead forcast error fce=abs(x(nw+1)-fcy) c Normalize x and => v call stat(lw,x(nw),am,sig,sk,c4,c6,smax,smin,io) rng=smax-smin if(sig>0.0) then do j=1,lw v(j)=(x(nn+j)-am)/sig enddo endif if(sig<=0.0) then flg=0.0 endif endif c------- c Farley-Hinich tests for mean & variance shifts call fhtest(v,lw,c4,prm,prv,io) if(flg==1.0) then c T23 test statistics if frame has var>0 call c234(v,lw,ng3,ng4,t4,h3,c2,nbc,ntc,io) if(q>0) then npa=q else npa=lw/4 endif h(kp)=h3 c(kp)=c2 h4(kp)=t4 a(kp)=am/gm s(kp)=sig/gs k3(kp)=sk k4(kp)=c4/g4 rg(kp)=rng/gr r2(kp)=ar2 c Check for variance & scale nonstationarity if(abs(s(kp)-out(6))>2.33*out(4)) then sht(kp)=s(kp) else sht(kp)=0. endif if(abs(a(kp)-out(6))>2.33*out(5)) then rht(kp)=a(kp) else rht(kp)=0. endif endif if(flg==0.0) then h(kp)=0. c(kp)=0. h4(kp)=0. a(kp)=0. s(kp)=0. k3(kp)=0. k4(kp)=0. rg(kp)=0. r2(kp)=0. sht(kp)=0. rht(kp)=0. endif c------- np2=np/2 c Test H statistic if(1.-h30.and.ov==0) then do j=1,lw w(j,nch)=x(nn+j) enddo endif else c Sum forecast for insignificant frame sumfn=sumfn+fce sefn=sefn+fce*fce endif c------- if(gf(1:1)=='y') then c Write statistics to graph file 4 if(q>0) then c AR fit c Check to see if time stamp read in. If mk='m' write it out. if(mk=='m') then write(4,'(sp,i7,1x,f7.4,12(1x,f11.4),1x,a20)')kp,h3,c2,t4,s(kp *),sht(kp),rg(kp),a(kp),rht(kp),k3(kp),k4(kp),r2(kp),prm,prv,t(nw) else write(4,'(sp,i7,1x,f7.4,12(1x,f11.4))')kp,h3,c2,t4,s(kp),sht( *kp),rg(kp),a(kp),rht(kp),k3(kp),k4(kp),r2(kp),prm,prv endif endif if(q==0) then c No AR if(mk=='m') then write(4,'(sp,i7,1x,f7.4,11(1x,f11.4),1x,a20)')kp,h3,c2,t4,s(kp *),sht(kp),rg(kp),a(kp),rht(kp),k3(kp),k4(kp),prm,prv,t(nw) else write(4,'(sp,i7,1x,f7.4,11(1x,f11.4))')kp,h3,c2,t4,s(kp),sh *t(kp),rg(kp),a(kp),rht(kp),k3(kp),k4(kp),prm,prv endif endif endif c End of writes to graph file enddo c End of loops on frames call flush(7) c------- if(gf(1:1)=='y') then if(q>0) then write(4,'(/2x,''Wind'',4x,''H'',11x,''C'',10x,''T4'',10x,''Sig'' *,9x,''Sig-CI'',6x,''Range'',7x,''Mean'',8x,''Mean-CI'',6x,''Skew'' *,8x,''Kurt'',7x,''R**2'',8x,''Mshift'',6x,''Vshift'')') else write(4,'(/2x,''Wind'',4x,''H'',11x,''C'',10x,''T4'',10x,''Sig *'',9x,''Sig-CI'',6x,''Range'',7x,''Mean'',8x,''Mean-CI'',6x,''Skew *'',8x,''Kurt'',7x,''Mshift'',6x,''Vshift'')') endif write(4,'(/2x,''Sig-CI are standard deviations that are outside a * 99% bootstrapped confidence interval'',/'' of its bootstrapped me *an'')') write(4,'(2x,''Ran-CI are ranges that are outside a 99% bootstrap *ped confidence interval of its'',/'' bootstrapped mean'')') write(4,'('' The mean (Mean), standard dev (Sig), kurtosis (Kurt) * and range for each frame are'')') write(4,'(''divided by their values computed from all the data to * normalize these statistics for plotting'')') write(4,'(4x,''Mshift is the Farley-Hinich test of a mean shift'' *)') write(4,'(4x,''Vshift is the Farley-Hinich test of a variance shi *ft'')') write(4,'(7x,''The values are cdf probabilities under the null of * no shift''/)') endif call flush(4) c------- if(gf(2:2)=='w'.and.ov==0) then c Output significant frames do j=1,nch write(11,'(sp,7000(:,f10.3,1x))')(w(k,j),k=1,lw) enddo endif c------- if(gf(1:1)=='y') then write(7,'(/9x,i4,'' Significant H frames'')')nch write(7,'(9x,i4,'' Significant C frames'')')ncc write(7,'(9x,i4,'' Significant T4 frames'')')nt4 endif if(nch>0) then write(io,'(/7x,i5,2x,''('',f6.2,''%)'','' significant H frames in * run'',i3/)')nch,100.*float(nch)/float(np),nr endif if(ncc>0) then write(io,'(7x,i5,2x,''('',f6.2,''%)'','' significant C frames''/) *')ncc,100.*float(ncc)/float(np) endif if(nt4>0) then write(io,'(7x,i5,2x,''('',f6.2,''%)'','' significant T4 frames''/ *)')nt4,100.*float(nt4)/float(np) endif if(ntt==0) then write(io,'(7x,''No signficant C, H or T4 statistics in run '',i4 *)')nr return endif c------ c Calculate correlations between stats if np>7 if(np<=7) then return endif c------ return end c c=========================================== c subroutine show(h3,c2,t4,c,e,r2,se,q,io) c************************************************* c Outputs statistics for data for a frame with significant stats c c - AR(q) coeff's, h3 - bicov stat, c2 - covariance stat, t4 - tricov stat c================================================= integer q real e(q),r2,se real*8 c(q) intent(in) h3,c2,t4,c,e,r2,se,q,io c************************************************* c Write out test statistics write(7,'(/4x,''p-values for tests: C = '',f9.7,2x,''H = '',f9.7,2 *x,''T4 = '',f9.7)') 1.-c2,1.-h3,1.-t4 if(q>0) then call disar(c,e,r2,se,q,7) endif c------- return end c c=========================================== c subroutine winhead(t,sr,lw,nn,unit,su,kn,mk,io) c******************************************* c Frame header for frames of length lw c kn - frame number, mk - switch for time mark, nn - pointer for frame start c=========================================== character*25 t(*) character*7 unit character su,mk intent(in) t,sr,lw,nn,unit,su,kn,mk,io c******************************************** rs=float(nn)+1. re=float(nn)+float(lw) write(io,'(70(''*''))') if(mk=='m') then write(io,'('' Frame '',i4,2x,''Start = '',a20,2x,'' End = '',a20 *)') kn,t(nn+1),t(nn+lw) elseif(su=='t') then write(io,'('' Frame '',i4,2x,''Start ='',f11.4,1x,a7,2x,'' End * ='',f11.4,1x,a7)') kn,rs,unit,re,unit elseif(su=='f') then del=1.e-3/sr write(io,'('' Frame '',i4,2x,''Start ='',g14.4,1x,'' sec'',2x, *''End ='',g14.4,1x,''sec'')') kn,rs*del,re*del else write(io,*) 'Flag mk for time file read is wrong ',mk return endif write(io,'(70(''=''))') c------- return end c c============================================ c subroutine cors(h,c,h4,s,k3,k4,rg,r2,indx,sumfs,sumfn,sefs,sefn,q, *np,lw,nch,ncc,io) c******************************************** c Correlations of test statistics c nch - no of frames with significant H statistics c sumfs,sefs - sum & sum of squares of absolute errors of forecasts for c signficicant H frames. sumfn,sefn - sum & sum of squares of absolute errors c of forecasts for insignficicant H frames. c============================================================ real h(np),c(np),h4(np),s(np),k3(np),k4(np),rg(np),r2(np) integer indx(np*lw),q,lw,np,nch,ncc,io character f1x,f1y,f2y intent(in) h,c,h4,s,k3,k4,rg,r2,sumfs,sefs,sumfn,sefn,q,np,lw,nch, *ncc,io c************************ smnfs=sumfs/float(nch) sigs=sefs/float(nch)-smnfs*smnfs if(sigs>0.0) then sigs=sqrt(sigs) endif c Trap nch if(nch==np) then sign=0.0 else smnfn=sumfn/float(np-nch) sign=sefn/float(np-nch)-smnfn*smnfn endif if(sign>0.0) then sign=sqrt(sign) endif c------ write(io,'('' Mean & Stnd. Deviation of One Step Forecast Absolute * Error for Significant Frames''/,2(f10.4,2x))')smnfs,sigs if(sign>0.0) then write(io,'(/'' Mean & Stnd Deviation of One Step Forecast Absolut *e Error for Insignificant Frames''/,2(f10.4,2x))')smnfn,sign endif write(io,'(2x,57(''*'')/)') c------- ish=np-nch+1 isc=np-ncc+1 ieh=nch-1 iec=ncc-1 c Correlations between significant H stats if nch>9 if(nch>9) then write(io,'(10x,''Correlations for '',i7,'' largest H Statistics'' *)')nch write(io,'(10x,42(''='')/)') call indexx(np,h,indx) call cor(np,ish,np,h,c,indx,rhc,f1x,f1y,io) call cor(np,ish,np,h,s,indx,rhs,f1x,f2y,io) c------- if(f1x=='y' .and. f1y=='y' .and. f2y=='y') then write(io,'(5x,''Correlation (H,C) = '',f5.2,7x,''Correlation (H,Si *g) = '',f5.2)')rhc,rhs endif if(f1x=='n') then write(io,*)'Var(H) = 0 for H frames' endif if(f1y=='n') then write(io,*)'Var(C) = 0 for H frames' endif if(f2y=='n') then write(io,*)'Var(Sig) = 0 for H frames' endif c------- if(q>0) then call cor(np,ish,np,h,r2,indx,hr2,f1x,f1y,io) call cor(np,ish,np,c,r2,indx,cr2,f1x,f2y,io) if(f1x=='y' .and. f1y=='y' .and. f2y=='y') then write(io,'(5x,''Correlation (H,R*R) = '',f5.2,7x,''Correlation (C, *R*R) = '',f5.2)')hr2,cr2 endif if(f1y=='n') then write(io,*)'Var(R*R) = 0 for H frames' endif endif c------- call cor(np,ish,np,h,k4,indx,rh4,f1x,f1y,io) call cor(np,ish,np,h,rg,indx,rhr,f1x,f2y,io) if(f1x=='y' .and. f1y=='y' .and. f2y=='y') then write(io,'(5x,''Correlation (H,Kurt) = '',f5.2,7x,''Correlation *(H,Range) = '',f5.2/)')rh4,rhr endif if(f1y=='n') then write(io,*)'Var(Kurt) = 0 for H frames' endif if(f2y=='n') then write(io,*)'Var(Range) = 0 for H frames' endif c------- write(io,'(10x,''Correlations of Statistics for the Rest of the W *indows'')') write(io,'(10x,54(''='')/)') call cor(np,1,ieh,h,c,indx,rhc,f1x,f1y,io) call cor(np,1,ieh,h,s,indx,rhs,f1x,f2y,io) if(f1x=='y' .and. f1y=='y' .and. f2y=='y') then write(io,'(5x,''Correlation (H,C) = '',f5.2,7x,''Correlation (H, *Sig) = '',f5.2)')rhc,rhs endif if(f1x=='n') then write(io,*)'Var(H) = 0 for the rest of the frames' endif if(f1y=='n') then write(io,*)'Var(C) = 0 for the rest of the frames' endif if(f2y=='n') then write(io,*)'Var(Sig) = 0 for the rest of the frames' endif c------- if(q>0) then call cor(np,ieh,np,h,r2,indx,hr2,f1x,f1y,io) call cor(np,ieh,np,c,r2,indx,cr2,f1x,f2y,io) if(f1x=='y' .and. f1y=='y' .and. f2y=='y') then write(io,'(5x,''Correlation (H,R*R) = '',f5.2,7x,''Correlation (C, *R*R) = '',f5.2)')hr2,cr2 endif if(f1y=='n') then write(io,*)'Var(R*R) = 0 for H frames' endif endif c------- call cor(np,1,ieh,h,k4,indx,rh4,f1x,f1y,io) call cor(np,1,ieh,h,rg,indx,rhr,f1x,f2y,io) if(f1x=='y' .and. f1y=='y' .and. f2y=='y') then write(io,'(5x,''Correlation (H,Kurt) = '',f5.2,7x,''Correlation *(H,Range) = '',f5.2/)')rh4,rhr endif if(f1y=='n') then write(io,*)'Var(Kurt) = 0 for H frames' endif endif c------- c Correlations between significant C stats if ncc>9 if(ncc>9) then write(io,'(10x,''Correlations for '',i4,'' largest C Statistics'' *)') ncc write(io,'(10x,42(''='')/)') c------- call indexx(np,c,indx) call cor(np,isc,np,c,h,indx,rhc,f1x,f1y,io) call cor(np,isc,np,c,s,indx,rcs,f1x,f2y,io) if(f1x=='y' .and. f1y=='y' .and. f2y=='y') then write(io,'(5x,''Correlation (C,H) = '',f5.2,7x,''Correlation (C, *Sig) = '',f5.2)')rhc,rcs endif if(f1x=='n') then write(io,*)'Var(C) = 0 for C frames' endif if(f1y=='n') then write(io,*)'Var(H) = 0 for C frames' endif if(f2y=='n') then write(io,*)'Var(Sig) = 0 for C frames' endif c------- if(q>0) then call cor(np,isc,np,c,r2,indx,cr2,f1x,f2y,io) if(f1x=='y' .and. f1y=='y' .and. f2y=='y') then write(io,'(5x,''Correlation (C,R**2) = '',f5.2)')cr2 endif if(f1y=='n') then write(io,*)'Var(R*R) = 0 for C frames' endif endif c------- call cor(np,isc,np,c,k4,indx,rc4,f1x,f1y,io) call cor(np,isc,np,c,rg,indx,rcr,f1x,f2y,io) if(f1x=='y' .and. f1y=='y' .and. f2y=='y') then write(io,'(5x,''Correlation (C,Kurt) = '',f5.2,7x,''Correlation *(C,Range) = '',f5.2/)')rc4,rcr endif if(f1y=='n') then write(io,*)'Var(Kurt) = 0 for C frames' endif if(f2y=='n') then write(io,*)'Var(Range) = 0 for C frames' endif c------- write(io,'(10x,''Correlations of Statistics for the Rest of the W *indows'')') write(io,'(10x,54(''='')/)') call cor(np,1,iec,c,h,indx,rhc,f1x,f1y,io) call cor(np,1,iec,c,s,indx,rhs,f1x,f2y,io) if(f1x=='y' .and. f1y=='y' .and. f2y=='y') then write(io,'(5x,''Correlation (C,H) = '',f5.2,7x,''Correlation (C, *Sig) = '',f5.2)')rhc,rhs endif if(f1x=='n') then write(io,*)'Var(C) = 0 for the rest of the frames' endif if(f1y=='n') then write(io,*)'Var(H) = 0 for the rest of the frames' endif if(f2y=='n') then write(io,*)'Var(Sig) = 0 for all frames' endif c------- if(q>0) then call cor(np,iec,np,c,r2,indx,cr2,f1x,f2y,io) if(f1x=='y' .and. f1y=='y' .and. f2y=='y') then write(io,'(5x,''Correlation (C,R**2) = '',f5.2/)')cr2 endif if(f1y=='n') then write(io,*)'Var(R*R) = 0 for the rest of the frames' endif endif c------- call cor(np,1,iec,c,k4,indx,rh4,f1x,f1y,io) call cor(np,1,iec,c,rg,indx,rhr,f1x,f2y,io) if(f1x=='y' .and. f1y=='y' .and. f2y=='y') then write(io,'(5x,''Correlation (C,Kurt) = '',f5.2,7x,''Correlation *(C,Range) = '',f5.2/)')rh4,rhr endif if(f1y=='n') then write(io,*)'Var(Kurt) = 0 for the rest of the frames' endif endif c------- return end c c============================================ c subroutine c234(x,lw,ng3,ng4,h4,h3,c,nbc,ntc,io) c******************************************** c Correlations, bicorrelations & tricorrelations test statistics for a frame c h4 - tranformed sum4 to uniform under the null (chi**2 with df=ntc) c h3 - tranformed sum3 to uniform under the null (chi**2 with df=nbc) c c - transformed sum2 to uniform under the null (chi**2 with df=ng3) c ng3 - no of lags for c & h3, ng4 - no of lags for h4 Version 7-21-2007 c============================================================ real x(lw) real*8 cv,bc,tc,sum2,sum3,sum4,nk1 integer t intent(in) x,lw,ng3,ng4,io intent(out) h4,h3,c,nbc,ntc c************************ c Trap lags if(ng4<3) then ntc=0 h4=1. return endif if(ng3<=ng4) then write(6,'(/''ERROR - Tricor lag '',i5,'' >= Bicor lag '',i5)')ng3 *,ng4 return endif c Initialize counters nbc=0 ntc=0 sum3=0.d0 sum4=0.d0 c------- c Compute cov(1) & b(1,2) cv=0.d0 bc=0.d0 do t=1,lw-2 cv=cv+dble(x(t))*dble(x(t+1)) bc=bc+dble(x(t))*dble(x(t+1))*dble(x(t+2)) enddo cv=cv+dble(x(lw-1))*dble(x(lw)) sum2=cv/dfloat(lw-1)*cv c------- c Start loop on k1 do k1=2,ng3 nk1=dfloat(lw-k1) c Loop on k20) then write(io,'(2x,''Statistics are computed from an AR('',i2,'') fit *to each frame'')')q endif c------- if(ov>0) then write(io,'(/2x,''Frames are overlapped by '',i6,'' observations'' *)')ov write(io,'(/'' Frame Length'',i5,2x,''No. of frames ='',i8)')lw *,np else write(io,'(/2x,''Frames are not overlapped'')') write(io,'(/'' Frame Length'',i5,2x,''No. of frames ='',i8)') *lw,np endif c------- if(cl=='tr'.or.cl=='tg') then write(io,'(/'' Lower & upper trim quantiles = '',2(f12.3,2x))')pc *l,pcu write(io,'(/'' Lower & upper trim levels = '',2(g12.3,2x),/3x,''N *umber Trimmed from Top & Bottom = '',2(i8,2x)/)')wl,wu,nuc,nlc write(io,'(60(''*''))') endif if(cl=='cl') then write(io,'(/7x,''Data is hard clipped around the mean'')') endif if(cl=='rr'.or.cl=='tg') then write(io,'(/7x,''Data transformed to growth rates)'')') endif write(io,'(/'' First datum read in file'',i8,4x,''Last datum read' *',i8)') is,ie write(io,'(/'' Sample size = '',i8/)') np*lw if(su=='t') then write(io,'('' Sampling interval = '',f12.6,2x,a7/)')sr,unit elseif(su=='f') then write(io,'('' Sampling rate = '',f10.4,'' kHz''/)')sr else write(io,*) 'Sample unit char in data file is not t or f' endif write(io,'(7x,''False alarm threshold = '',f7.5/)')th write(io,'(4x,''C & H Exponent = '',f5.2,2x,''T Exponent = '',f5.2 */)')exp,te ng3=nint(float(lw)**exp) ng4=nint(float(lw)**te) write(io,'(4x,''C & H Lags = '',i6,2x,''T4 Lags = '',i6/)')ng3,ng4 write(io,'(4x,''Bicovariances used = '',i9,2x,''Tricoviances used *= '',i8/)')nbc,ntc if(res>0) then write(io,'(7x,i4,'' Bootstrapped Sizes: H = '',g9.3,2x,''C = '',g *9.3,2x,''T4 = '',g9.3/)')res,out(1),out(2),out(3) write(io,'(11x,''Bootstrapped Standard Errors''/)') write(io,'(7x,''For Relative Standard Deviations = '',f10.4/)')ou *t(4) write(io,'(7x,''For Relative Ranges = '',f10.4/)')out(5) write(io,'(7x,''For Relative Kurtoses = '',f10.4/)')out(6) endif write(io,'('' Format: '',a50)')fmt return end c c=================================================== c subroutine wrstat(i,am,sig,sk,c4,c6,smax,smin) c**************************************************** write(i,'(60(''='')/)') write(i,'('' Mean ='',g12.3,7x,''Std Dev ='',g12.3/)') am,sig write(i,'('' Skew ='',g10.3,4x,''Kurtosis ='',g10.3,4x,''C(6) ='', *f10.2/)') sk,c4,c6 write(i,'('' Max value ='',g14.3,7x,''Min value ='',g14.3)') smax, *smin write(i,'(60(''*'')/)') return end c c=========================================== c subroutine datc(x,t,fmt,mk,is,ie,n,nc,ii,io) c******************************************************** c Reads n data values from nc columns c x - input data matrix, t - time marks if indicated by 'm' c======================================================== real x(n,nc) character*50 fmt character*25 t(*) character mk intent(in) fmt,mk,is,ie,n,nc,ii,io intent(out) x,t c******************************************** nt=0 if(mk=='m') then ia=index(fmt,'(') ib=index(fmt,'a') if(ib==ia+1) then c Time mark column on left nt=1 elseif(ib>ia+2) then c Time mark column on right nt=2 elseif(nt==0) then write(io,'(7x,''/Data & time format is incorrect''/)') write(io,'(a50)')fmt write(6,'(7x,''/Data & time format is incorrect''/)') write(6,'(a50)')fmt write(9,'(7x,''/Data & time format is incorrect''/)') write(9,'(a50)')fmt stop endif endif c------- rewind(ii) c Skip header and parameter records read(ii,*) read(ii,*) c Read all the data if(mk=='m') then if(is==1) then do i=1,ie if(nt==1) then read(ii,fmt,end=1,err=2,iostat=ko)t(i),(x(i,k),k=1,nc) endif if(nt==2) then read(ii,fmt,end=1,err=2,iostat=ko)(x(i,k),k=1,nc),t(i) endif enddo endif if(is>1) then c Read is:ie do i=1,is-1 read(ii,fmt,end=1,err=2,iostat=ko) enddo do i=is,ie if(nt==1) then read(ii,fmt,end=1,err=2,iostat=ko)t(i-is+1),(x(i-is+1,k),k=1,n *c) endif if(nt==2) then read(ii,fmt,end=1,err=2,iostat=ko)(x(i-is+1,k),k=1,nc),t(i-is+ *1) endif enddo endif endif if(mk/='m') then if(is==1) then do i=1,ie read(ii,fmt,end=1,err=2,iostat=ko)(x(i,k),k=1,nc) enddo endif if(is>1) then do i=1,is-1 read(ii,fmt,end=1,err=2,iostat=ko) enddo do i=is,ie read(ii,fmt,end=1,err=2,iostat=ko)(x(i-is+1,k),k=1,nc) enddo endif endif return c------- 1 write(io,'(/2x,''ERROR End of file in reading the data file'')') write(io,'(/7x,'' at data index t = '',i6)')i write(6,'(/2x,''ERROR End of file in reading the data file'')') write(6,'(/7x,'' at data index t = '',i6)')i write(9,'(/2x,''ERROR End of file in reading the data file'')') write(9,'(/7x,'' at data index t = '',i6)')i stop c Error in data read 2 write(io,*) '** Error in read of data file ',ko write(io,*) 'Check format',fmt write(6,*) '** Error in read of data file ',ko write(6,*) 'Check format',fmt write(9,*) '** Error in read of data file ',ko write(9,*) 'Check format',fmt stop end c c============================================ c subroutine forecast(y,c,ymn,fcy,n,p) c******************************************** c One step ahead forecast of AR(p) c============================================ integer p,t real y(n),ymn real*8 c(p),sm intent(in) y,c,ymn,n,p intent(out) fcy c******************************************** sm=0.d0 do t=1,p sm=sm+c(t)*(dble(y(n+1-t)-ymn)) enddo fcy=sngl(sm)+ymn c------- return end c c=========================================== c subroutine fhtest(x,n,c4,prm,prv,io) c******************************************** c Farley-Hinich tests for mean & variance shifts c prm - p-value for mean shift, prv - p-value for variance shift c============================================ real x(n) real*8 dvn,sumx,sumv,vt,xt,zn integer t intent(in) x,n,c4,io intent(out) prm,prv c**************************** if(dble(c4)+2.d0<=0.d0) then prm=1.0 prv=1.0 return endif c Standard deviation of sum((t-(n+1)/2)*x(t) hn=float(n+1)/2. s2n=float(2*n+1)/3. c sqrt(n*(n+1)*(2n+1)/6) sdtsum=sqrt(float(n)*hn*s2n) dvn=dble(hn) c Mean test sumx=0.d0 sumv=0.d0 c Standard deviation of x(t)**2 zn=dsqrt(dble(c4)+2.d0) do t=1,n c Centered t tn=dfloat(t)-dvn xt=dble(x(t)) c Standardized x(t)**2 vt=(xt*xt-1.d0)/zn sumx=sumx+tn*xt sumv=sumv+tn*vt enddo c P-values prm=1.-cdg(sngl(sumx)/sdtsum) prv=1.-cdg(sngl(sumv)/sdtsum) c------- return end c c=========================================== c subroutine resamp(x,y,n,lw,res,exp,te,gm,gs,g4,thrsh,out,ir,io) c******************************************** c 1% & 5% bootstrap thresholds & standard deviations for means, sigs, & c kurtosis for resamples=res, ir=-1 - 0 sig y Version 4-23-2008 c============================================ allocatable::ah,ac,ah4,asg,arg,ac4,indx real x(n),y(lw),ah(:),ac(:),ah4(:),asg(:),arg(:),ac4(:),fi(1),fo(1 *),out(6) integer indx(:),res,t intent(in) x,n,lw,res,exp,te,gm,gs,g4,thrsh,io intent(out) out,ir c******************************************** c Allocate space ns=res-1 nw=lw-1 allocate(ah(0:ns),ac(0:ns),ah4(0:ns),asg(0:ns),arg(0:ns),ac4(0:ns) *,indx(0:ns),stat=ibad) if(ibad/=0) then write(io,*) 'Unable to allocate work space for arrays in Resamp' write(6,*) 'Unable to allocate work space for arrays in Resamp' stop endif c------ rn=float(n) do kr=0,res-1 call random_number(y) do t=1,lw nsam=max(nint(y(t)*rn),1) y(t)=x(nsam) enddo ir=0 call stat(lw,y,am,sig,sk,c4,c6,smax,smin,io) if(sig<1.e-9) then ir=-1 return endif do t=1,lw y(t)=(y(t)-am)/sig enddo c Normalized mean arg(kr)=am/gm c Normalized sigima asg(kr)=sig/gs c Normalized c4 ac4(kr)=c4/g4 c C234 of y ng3=nint(float(lw)**exp) ng4=nint(float(lw)**te) call c234(y,lw,ng3,ng4,t4,h3,c2,nbc,ntc,io) ah(kr)=1.0-h3 ac(kr)=1.0-c2 ah4(kr)=1.0-t4 enddo c End of resample loop c------- c Size threshold fi(1)=thrsh call sort(indx,ah,res,1,fi,fo,'B',ier) c Threshold for qv% h3, c2 & t4 out(1)=fo(1) call sort(indx,ac,res,1,fi,fo,'B',ier) out(2)=fo(1) call sort(indx,ah4,res,1,fi,fo,'B',ier) out(3)=fo(1) c Standard deviations of sigmas, means, c4s call sigmn(asg,res,am,sigs,io) out(4)=sigs call sigmn(arg,res,am,sigm,io) if(sigm>0.0) then out(5)=sigm endif if(sigm==0.0) then out(5)=thrsh endif call sigmn(ac4,res,am,sig4,io) out(6)=sig4 c------- deallocate(ah,ac,ah4,asg,arg,ac4,indx) return end c c=================================================== c subroutine detrend(y,r,c,in,n,pv,po,ram,rsd,io) c******************************************************** c Detrend using an ols fit of 1-4 th order Legendre polynomials 11-4-2007 c pv - p Value Threshold for the Iterative OLS fit c======================================================== implicit none allocatable::x real x(:,:),y(n),r(n),pv,rn,z,z2,z3,z4,r2,ram,rsd,seb,rsk,rc4,rc6, *sig,ts,am,rmax,rmin real*8 a(16),c(4) integer in(4),po,t,n,n2,k,ibad,ier,io character*2 varnm(4) intent(in) y,n,pv,io intent(out) r,c,in,po,ram,rsd c********************** c Allocate space allocate(x(0:n,0:3),stat=ibad) if(ibad/=0) then write(io,*)'Unable to allocate work space for array in detrend' write(6,*)'Unable to allocate work space for array in detrend' return endif c------- c Legendre polynomials n2=n/2 rn=float(n2) do t=-n2,n2 z=float(t)/rn z2=z*z z3=z*z2 z4=z*z3 k=t+n2 x(k,0)=z x(k,1)=(3.0*z2-1.0)/2.0 x(k,2)=(5.0*z3-3.0*z)/2.0 x(k,3)=(35.0*z4-30.0*z2+3.0)/8.0 enddo varnm(1)='L1' varnm(2)='L2' varnm(3)='L3' varnm(4)='L4' c Least squares fit of y=x*b+e ier=0 call ols(y,x,c,r,a,in,am,r2,sig,n,4,po,pv,ier,io) if(ier>0) then write(io,*)'Error in subroutine ols ier = ',ier write(6,*)'Error in subroutine ols ier = ',ier po=0 ram=0.0 rsd=1.0 return endif if(po==0) then write(io,'(/4x,''The Trend is Statistically Insignificant'')') write(io,'(7x,''p Value Threshold for the Iterative OLS Fit = '', *f5.3)')pv ram=0.0 rsd=1.0 return endif c------- write(io,'(/7x,''Data is Detrended by a Least Squares Fit of Four *Legendre Polynomials'')') write(io,'(7x,''Significant Independent Variables = '',i3/)')po write(io,'(7x,''p Value Threshold for the Iterative OLS Fit = '',f *5.3)')pv write(io,'(1x,64(''-''))') write(io,'(/12x,''Coefficients & Estimates'')') write(io,'(1x,64(''*''))') c Write out estimates do k=1,po seb=sig*sngl(dsqrt(a((k-1)*po+k))) ts=sngl(c(k))/seb write(io,'(2x,''b^('',a2,'') = '',f7.3,2x,''Std Error = '',f9.5,2 *x,''t = '',f10.3)')varnm(in(k)),c(k),seb,ts enddo write(io,'(64(''=''))') c Write R square write(io,'(/3x,''Adjusted R Squared = '',f7.4,2x,''Std Error of Re *siduals= '',f9.4)')r2,sig write(io,'(/7x,''Estimate of Intercept = '',f9.3/)')am c Statistics call stat(n,r,ram,rsd,rsk,rc4,rc6,rmax,rmin,io) write(io,'(57(''=''))') write(io,'(7x,''Statistics for the Residuals'')') write(io,'(57(''*'')/)') write(io, '(7x,''Mean = '',g12.4,2x,''Sig = '',g12.4/)')ram,rsd write(io,'(7x,'' Skewness = '',g12.3,2x,''Kurtosis = '',f12.2/)') *rsk,rc4 write(io,'(7x,'' Max ='',g17.5,2x,''Min ='',g17.5)')rmax,rmin c------- deallocate(x) return end c c=================================================== c subroutine trim(x,y,n,pcu,pcl,wu,wl,io) c**************************************************** c Melvin Hinich Version 11-2-2006 c Trim data in array x at pcl% & pcu% quantiles c================================== c Input: x - input data array, pcl, pcu - % trimmed from extremes of df c Output: y - clipped data array, wu & wi - top & bottom pc/2 quantiles c================================== allocatable::indx real x(n),y(n),frac(1),vfrac(1) integer indx(:),t intent(in) x,n,pcu,pcl intent(out) y c********************************** if(pcu>=100.0) then write(6,*)'Upper % trim value out of bounds in Trim ',pcu return endif if(pcl<=0.0) then write(6,*)'Lower % trim value out of bounds in Trim ',pcl return endif c------- allocate(indx(0:n-1),stat=ibad) if(ibad/=0) then write(io,*)'Unable to allocate work space for indx in trim' write(6,*)'Unable to allocate work space for indx in trim' stop endif c------- c Compute number to be changed at each end nlc=nint(pcl/100.*float(n)) nuc=nint((1.0-pcu/100.0)*float(n)) if(nuc==0) then nuc=1 endif if(nlc==0) then nlc=1 endif do t=1,n y(t)=x(t) enddo c Sort the wk array call sort(indx,x,n,1,frac,vfrac,'S',ier) if(ier>0) then write(6,'(/7x,''Error in Sort in Trim''/)') stop endif c Datum at the upper quantile wu=x(indx(n-nuc-1)) if(nlc>0) then c Datum at the lower quantile wl=x(indx(nlc-1)) else wl=x(indx(0)) endif range=wu-wl if(range<=0.0) then write(io,'(/7x,''Range '',f10.2,'' in subroutine Trim <= 0''/)')r *ange return endif c------- do t=1,n if(x(t)>wu) then c Trim large y(t)=wu elseif(x(t)y y(t)=x(t) endif enddo c------- deallocate(indx) return end c c=================================================== c subroutine areturns(y,r,n,io) c******************************************** c Rates of returns => r c============================================ allocatable::x2 real y(n),r(n),x2(:),yy integer t,n1 intent(in) y,n,io intent(out) r c******************************************* c Allocate space allocate(x2(0:n-1),stat=ibad) if(ibad/=0) then write(io,*)'Unable to allocate work space for arrays in returns' write(6,*)'Unable to allocate work space for arrays in returns' stop endif c------- c Copy source data to scratch array x2 do t=1,n x2(t-1)=0.0 enddo x2(0)=y(1) x2(n-1)=y(n) do t=1,n-2 x2(t)=y(t+1) enddo c Calculate returns r(1)=0.0 do t=1,n-1 yy=x2(t-1) if(yy<=0.0) then r(t+1)=0.0 elseif((abs(x2(t)/yy))<1.e-17) then write(io,'(2x,''return in Returns < 1.e-17 for t = '',i6)')t write(io,'(2x,''return set to zero'')') r(t+1)=0.0 else r(t+1)=100.0*log(x2(t)/yy) endif enddo deallocate(x2) c------- return end c c=================================================== c subroutine hardclp(y,r,n,mean,io) c******************************************** c Hard clipped data => r c============================================ implicit none real y(n),r(n),mean integer n,t,io intent(in) y,n,mean intent(out) r c******************************************** do t=1,n if(y(t)>=mean) then r(t)=1.0 else r(t)=-1.0 endif enddo c------- return end