C00000000111111111122222222223333333333444444444455555555556666666666777
C23456789012345678901234567890123456789012345678901234567890123456789012
C
C  THESE ROUTINES CALCULATE THE COVARIANCE MATRIX FOR THE PARAMETER
C  ESTIMATES
C***********************************************************************
C       COVAR---CALCULATE THE COVARIANCE MATRIX
C       JACOB---FIND THE JACOBIAN FOR MIXED ROOTS
C       SIGUV---TERM BY TERM CALCULATION FOR COVARIANCES IN THE FORM
C               FROM DR. LII'S PAPER  NOTE:  ERRORS FOUND WITH
C               TERM 8 AND SIGN CHANGES WITH TERMS 14 AND 15
C       SHIFT---A ROUTINE TO PLACE THE POLYNOMIALS IN SAME FORM AS
C               DESCRIBED IN DR. LII'S PAPER
C***********************************************************************
C
C       CALCULATE THE COVARIANCE MATRIX
C***********************************************************************
C       INPUT
C               C1,C2--THE COEFFICIENTS FROM THE DISTRIBUTION
C                      DESCRIBED IN DR. LII'S PAPER
C               COEFF--THE ARMA COEFFICIENTS
C               NCPHI--THE P IN ARMA(P,Q)
C               NCTHET-THE Q IN ARMA(P,Q)
C               MAXL----MAXIMUM LENGTH OF INVERSE POLYNOMIALS
C               ACURCY-THE ACCURACY OF THE INVERSE POLYNOMIALS
C               NSAMP--NUMBER OF SAMPLES
C               IQ-----TRUNCATION OF SERIES
C               NP-----DIMENSION OF POLYNOMIALS
C       OUTPUT
C               PARM---THE NON-ZERO COEFFICIENT PARAMETERS
C               COV----THE COVARIANCE MATRIX
C               STD----PARAMTER STANDARD DEVIATIONS
C               NC-----NUMBER OF NON-ZERO PARAMETERS
C***********************************************************************

      SUBROUTINE COVAR(C1,C2,COEFF,NCPHI,NCTHET,MAXL,ACURCY,
     1                 NSAMP,IQ,NP,PARM,COV,STD,NC,
     2                 NPI,NPO,NTI,NTO,PI,PO,TI,TO,
     3                 A,B,AP,BP,CV,RJ,RJT,EINFO,RLJ,RJCV,Y,PATH)
        INTEGER NCPHI,NCTHET,NP,INFO,MAXL,IQ,NC,ISIG,NETA,ICOUNT
        INTEGER I,JCOUNT,J,NSAMP
        INTEGER R,RP,S,SP,P,PP,IPATH
        INTEGER NA,NAP,NB,NBP,NPO,NPI,NTO,NTI
        REAL COEFF(NP),PARM(NP),COV(NP,NP),STD(NP)
        REAL A(0:MAXL),B(0:MAXL),AP(0:MAXL),BP(0:MAXL)
        REAL CV(NP,NP),RJ(NP,NP),RJT(NP,NP),EINFO(NP,NP)
        REAL RLJ(NP,NP),RJCV(NP,NP),PI(NP),PO(NP),TI(NP),TO(NP)
        REAL Y(NP,NP),PATH(NP,NP)
        REAL C1,C2,ACURCY,SIG,SIGUV
        CHARACTER UPLO

C  SET THE SCALE FACTOR
        ISIG=NCPHI+NCTHET+1
        SIG=COEFF(ISIG)

C  PUT INTO FORM FOR USE WITH DR.LII'S ASYMPTOTIC COVARIANCE MATRIX
        CALL SHIFT(PI,NPI,PO,NPO,TI,NTI,TO,NTO,NP,MAXL,ACURCY,
     1             R,RP,S,SP,P,PP,NETA,A,NA,AP,NAP,B,NB,BP,NBP)

C  FIND THE JACOBIAN
        CALL JACOB(PI,NPI,PO,NPO,TI,NTI,TO,NTO,NP,RLJ)

C  CALCULATE THE SIGMA MATRIX FROM DR. LII'S PAPER
C  ONLY CALCULATE THE NONZERO TERMS
        ICOUNT=0
        DO 20 I=1,NETA
                IF(COEFF(I).NE.0) THEN
                ICOUNT=ICOUNT+1
                PARM(ICOUNT)=COEFF(I)
                        JCOUNT=0
                        DO 10 J=1,NETA
                                IF(COEFF(J).NE.0) THEN
                                JCOUNT=JCOUNT+1
                RJ(ICOUNT,JCOUNT)=RLJ(I,J)
                EINFO(ICOUNT,JCOUNT)=SIGUV(C1,C2,I,J,R,RP,S,SP,P,PP,
     1                A,NA,AP,NAP,B,NB,BP,NBP,MAXL,SIG,IPATH)
                PATH(ICOUNT,JCOUNT)=IPATH
                ENDIF
10              CONTINUE
                ENDIF
20      CONTINUE

C  USE THE LAPACK ROUTINE TO INVERT THE EINFO MATRIX
        DO 40 I=1,ICOUNT
               DO 30 J=1,JCOUNT
                      Y(I,J)=0.0
30             CONTINUE
               Y(I,I)=1.0
40      CONTINUE
        UPLO='L'
        CALL SPOSV(UPLO,ICOUNT,ICOUNT,EINFO,NP,Y,NP,INFO)

C  CALCULATE THE INFORMATION MATRIX
        CALL MATMUL(RJ,ICOUNT,JCOUNT,Y,ICOUNT,JCOUNT,RJCV,NP)
        CALL TRNSPS(RJ,ICOUNT,JCOUNT,RJT,NP)
        CALL MATMUL(RJCV,ICOUNT,JCOUNT,RJT,ICOUNT,JCOUNT,CV,NP)

C  NOW CALCULATE THE COVARIANCE MATRIX AND THE PARAMETER STANDARD
C  DEVIATIONS

        DO 70 I=1,ICOUNT
        STD(I)=SQRT(CV(I,I)/(NSAMP-2*IQ))
               DO 60 J=1,JCOUNT
               COV(I,J)=CV(I,J)/(NSAMP-2*IQ)
60             CONTINUE
70      CONTINUE
        NC=ICOUNT

        RETURN
        END

C***********************************************************************
C       CALCULATE THE JACOBIAN USING THE POLYNOMIALS INSIDE AND
C       OUTSIDE THE UNIT MATRIX OF THE AR AND MA COEFFICIENTS
C***********************************************************************
C       INPUT
C               PI,PO----AR POLYNOMIALS INSIDE AND OUTSIDE UNIT CIRCLE
C               NPI,NPO--DEGREE OF ABOVE POLYNOMIALS
C               TI,TO----MA POLYNOMIALS INSIDE AND OUTSIDE UNIT CIRCLE
C               NTI,NTO--DEGREE OF ABOVE POLYNOMIALS
C               NP-------THE DIMENSION OF THE COEF POLYNOMIALS
C       OUTPUT
C               R--------THE JACOBIAN
C***********************************************************************

      SUBROUTINE JACOB(PI,NPI,PO,NPO,
     2                        TI,NTI,TO,NTO,NP,R)
        INTEGER NPI,NPO,NTI,NTO,NP,I,J,INDEX,INDEXR,NSIZE
        REAL PI(NP),PO(NP),TI(NP),TO(NP)
        REAL R(NP,NP)

C  INITILIZE THE JACOBIAN
        DO 20 I=1,NP
               DO 10 J=1,NP
                      R(I,J)=0
10             CONTINUE
20      CONTINUE

C  FILL THE JACOBIAN WITH THE PROPER VALUES  PHIOUT WHICH ARE THE
C  DERIVATIVES FOR PHIIN

        DO 40 I=1,NPO+NPI
               DO 30 J=1,NPO
                      INDEX=I-J+1
           IF((INDEX.GT.0).AND.(INDEX.LE.NPI+1)) R(I,J)=PI(INDEX)
30             CONTINUE
40      CONTINUE

C  PHIIN WHICH THE DERIVATIVES ARE PHIOUT
        DO 60 I=1,NPI+NPO
               DO 50 J=1,NPI
                      INDEX=I-J+1
        IF((INDEX.GT.0).AND.(INDEX.LE.NPO+1)) R(I,NPO+J)=PO(INDEX)
50             CONTINUE
60      CONTINUE

C  THETAOUT WHICH ARE THE DERIVATIVES OF THETAIN
        INDEXR=NPI+NPO
        DO 80 I=1,NTO+NTI
               DO 70 J=1,NTO
                      INDEX=I-J+1
        IF((INDEX.GT.0).AND.(INDEX.LE.NTI+1))
     1                R(INDEXR+I,INDEXR+J)=TI(INDEX)
70             CONTINUE
80      CONTINUE

C  THETAIN WHICH ARE THE DERIVATIVES GIVING THETAOUT
        DO 100 I=1,NTI+NTO
                DO 90 J=1,NTI
                       INDEX=I-J+1
                IF((INDEX.GT.0).AND.(INDEX.LE.NTO+1)) THEN
                       R(INDEXR+I,INDEXR+NTO+J)=TO(INDEX)
                ENDIF
90              CONTINUE
100     CONTINUE

C  THE JACOBIAN WITH RESPECT TO SIGMA
        NSIZE=NPI+NPO+NTI+NTO+1
        R(NSIZE,NSIZE)=1.0

        RETURN
        END

C***********************************************************************
C       CALCULATE THE INFORMATION MATRIX W.R.T. DR. LII'S PAPER
C       NOTE:  EXACTLY THE SAME TERMINOLOGY IS USED
C***********************************************************************
C       INPUT
C               C1,C2------DENSITY COEFFICIENTS
C               UI,VI------THE ROW AND COLUMN IN THE MATRIX
C       REFER TO SUBROUTINE SHIFT FOR DEFINITION OF THE OTHER
C       PARMETERS
C***********************************************************************

      FUNCTION SIGUV(C1,C2,UI,VI,R,RP,S,SP,P,PP,
     1                A,NA,AP,NAP,B,NB,BP,NBP,MAXL,SIG,IPATH)
        INTEGER UI,VI,U,V,R,RP,S,SP,P,PP,IPATH,MAXL,NA,NAP,NB,NBP
        INTEGER J,J1,J2,NL
        REAL A(0:MAXL),AP(0:MAXL),B(0:MAXL),BP(0:MAXL)
        REAL C1,C2,SIG,SUM,SIGUV

C  INITIALIZE PARAMETERS
        IPATH=0
        IF(UI.LE.VI) THEN
                U=UI
                V=VI
        ELSE
                U=VI
                V=UI
        ENDIF

C  SUM USED IN ALL TERMS BELOW
        SUM=0
        SIGUV=0


C  FIRST TERM
        IF((U.LE.R).AND.(V.LE.R)) THEN
                DO 30 J=0,NA
                        J1=J+ABS(U-V)
                        IF((J1.GE.0).AND.(J1.LE.NA)) THEN
                        SUM=SUM+A(J)*A(J1)                
                        ENDIF
30              CONTINUE
                SIGUV=C1*SUM
                IPATH=1
C  SECOND TERM
        ELSEIF((U.GT.R).AND.(U.LE.P).AND.
     1               (V.GT.R).AND.(V.LE.P).AND.
     1         (.NOT.((U.EQ.P).AND.(V.EQ.P)))) THEN
               DO 40 J=S,NB
                      J1=J+ABS(U-V)
                      IF((J1.GE.S).AND.(J1.LE.NB)) THEN
                        SUM=SUM+B(J)*B(J1)
                      ENDIF
40              CONTINUE
                SIGUV=C1*SUM
                IPATH=2
C  THIRD TERM
        ELSEIF((U.GT.R).AND.(U.LE.P).AND.
     1               (V.GT.R).AND.(V.LE.P).AND.
     1        ((U.EQ.P).AND.(V.EQ.P))) THEN
                DO 50 J=S,NB
                       IF((J+1).LE.NB) THEN
                       SUM=SUM+B(J+1)*B(J+1)  
                       ENDIF
50              CONTINUE
                SIGUV=B(S)*B(S)*(C2-1)+C1*SUM
                IPATH=3
C  FOURTH TERM
         ELSEIF((U.LE.R).AND.(V.GT.R).AND.(V.LE.P)) THEN
                IF(NA.GT.NB) THEN
                        NL=NA
                ELSE
                        NL=NB
                ENDIF
                DO 60 J=-NL,NL
                        J1=J-U
                        J2=J+V-R
                        IF((J1.GE.0).AND.(J1.LE.NA).AND.
     1                         (J2.GE.S).AND.(J2.LE.NB)) THEN
                               SUM=SUM+A(J1)*B(J2)
                        ENDIF
60              CONTINUE
                SIGUV=SUM
                IPATH=4
C  FIFTH TERM
        ELSEIF((U.GT.P).AND.(U.LE.(P+RP)).AND.
     1          (V.GT.P).AND.(V.LE.(P+RP))) THEN
                DO 70 J=0,NAP
                        J1=J+ABS(V-U)
                        IF((J1.GE.0).AND.(J1.LE.NAP)) THEN
                                SUM=SUM+AP(J)*AP(J1)
                        ENDIF
70              CONTINUE
                SIGUV=C1*SUM
                IPATH=5
C  SIXTH TERM
        ELSEIF((U.GE.(P+RP+1)).AND.(U.LE.(P+PP)).AND.
     1          (.NOT.((U.EQ.(P+PP)).AND.(V.EQ.(P+PP)))).AND.
     1          (V.GE.(P+RP+1)).AND.(V.LE.(P+PP))) THEN
                DO 80 J=SP,NBP
                        J1=J+ABS(V-U)
                        IF((J1.GE.SP).AND.(J1.LE.NBP)) THEN
                                SUM=SUM+BP(J)*BP(J1)
                        ENDIF
80              CONTINUE
                SIGUV=C1*SUM
                IPATH=6
C  SEVENTH TERM
        ELSEIF((U.GE.(P+RP+1)).AND.(U.LE.(P+PP)).AND.
     1          ((U.EQ.(P+PP)).AND.(V.EQ.(P+PP))).AND.
     1          (V.GE.(P+RP+1)).AND.(V.LE.(P+PP))) THEN
                DO 90 J=SP,NBP
                        IF(((J+1).GE.SP).AND.((J+1).LE.NBP)) THEN
                                SUM=SUM+BP(J+1)*BP(J+1)
                        ENDIF
90              CONTINUE
                SIGUV=BP(SP)*BP(SP)*(C2-1)+C1*SUM
                IPATH=7
C  EIGTH TERM  (NOTE:  THIS IS A CORRECTION TO THE PAPER)
        ELSEIF((U.GT.P).AND.(U.LE.(P+RP)).AND.
     1          (V.GT.(P+RP)).AND.(V.LE.(P+PP))) THEN
                IF(NAP.GT.NBP) THEN
                        NL=NAP
                ELSE
                        NL=NBP
                ENDIF
                DO 100 J=-NL,NL
                        J1=J-(U-P)
                        J2=J+V-P-RP
                        IF((J1.GE.0).AND.(J1.LE.NAP).AND.
     1                     (J2.GE.SP).AND.(J2.LE.NBP)) THEN
                                SUM=SUM+AP(J1)*BP(J2)
                        ENDIF
100             CONTINUE
                SIGUV=SUM 
                IPATH=8
C  NINTH TERM
        ELSEIF((U.LE.R).AND.(V.GT.P).AND.(V.LE.(P+RP))) THEN
                IF(NA.GT.NAP) THEN
                        NL=NA
                ELSE
                        NL=NAP
                ENDIF
                DO 110 J=-NL,NL
                        J1=J+U-V+P
                        IF((J.GE.0).AND.(J.LE.NA).AND.
     1                     (J1.GE.0).AND.(J1.LE.NAP)) THEN
                                  SUM=SUM+A(J)*AP(J1)
                        ENDIF
110             CONTINUE
                SIGUV=-C1*SUM
                IPATH=9
C  TENTH TERM
        ELSEIF((U.GT.R).AND.(U.LE.P).AND.
     1          (.NOT.((U.EQ.P).AND.(V.EQ.(P+PP)))).AND.
     1          (V.GT.(P+RP)).AND.(V.LE.(P+PP))) THEN
                IF(NB.GT.NBP) THEN
                        NL=NB
                ELSE
                        NL=NBP
                ENDIF
                DO 120 J=-NL,NL
                        J1=J+R-U-RP+V-P
                        IF((J.GE.S).AND.(J.LE.NB).AND.
     1                     (J1.GE.SP).AND.(J1.LE.NBP)) THEN
                                  SUM=SUM+B(J)*BP(J1)
                        ENDIF
120             CONTINUE
                SIGUV=-C1*SUM
                IPATH=10

C  ELEVENTH TERM
        ELSEIF((U.GT.R).AND.(U.LE.P).AND.
     1         ((U.EQ.P).AND.(V.EQ.(P+PP))).AND.
     1          (V.GT.(P+RP)).AND.(V.LE.(P+PP))) THEN
                IF(NB.GT.NBP) THEN
                        NL=NB
                ELSE
                        NL=NBP
                ENDIF
                DO 130 J=-NL,NL
                        J1=S+J
                        J2=SP+J
                        IF((J1.GT.S).AND.(J1.LE.NB).AND.
     1                     (J2.GT.SP).AND.(J2.LE.NBP)) THEN
                                 SUM=SUM+B(J1)*BP(J2)
                        ENDIF
130             CONTINUE
                SIGUV=-B(S)*BP(SP)*(C2-1)-C1*SUM
                IPATH=11
C  TWELVTH TERM
          ELSEIF((U.LE.R).AND.(V.GT.(P+RP)).AND.(V.LE.(P+PP))) THEN
                IF(NA.GT.NBP) THEN
                        NL=NA
                ELSE
                        NL=NBP
                ENDIF
                DO 140 J=-NL,NL
                        J1=J-U
                        J2=J+V-P-RP
                        IF((J1.GE.0).AND.(J1.LE.NA).AND.
     1                     (J2.GE.SP).AND.(J2.LE.NBP)) THEN 
                              SUM=SUM+A(J1)*BP(J2)
                        ENDIF
140             CONTINUE
                SIGUV=-SUM
                IPATH=12
C  THIRTEENTH TERM
        ELSEIF((U.GT.R).AND.(U.LE.P).AND.
     1         (V.GT.P).AND.(V.LE.(P+RP))) THEN
                IF(NAP.GT.NB) THEN
                        NL=NAP
                ELSE
                        NL=NB
                ENDIF
                DO 150 J=-NL,NL
                        J1=J-V+P
                        J2=J+U-R
                        IF((J1.GE.0).AND.(J1.LE.NAP).AND.
     1                     (J2.GE.S).AND.(J2.LE.NB)) THEN        
                                  SUM=SUM+AP(J1)*B(J2)
                        ENDIF        
150             CONTINUE
                SIGUV=-SUM
                IPATH=13 
C  FOURTEENTH TERM
        ELSEIF((U.EQ.P).AND.(V.EQ.(P+PP+1))) THEN
                SIGUV=-B(S)*(C2-1)/SIG
                IPATH=14
C  FIFTEENTH TERM
        ELSEIF((U.EQ.(P+PP)).AND.(V.EQ.(P+PP+1))) THEN
                SIGUV=+BP(SP)*(C2-1)/SIG
                IPATH=15
C  SIXTEENTH TERM
        ELSEIF((U.EQ.(P+PP+1)).AND.(V.EQ.(P+PP+1))) THEN
                SIGUV=(C2-1)/SIG/SIG
                IPATH=16
C SEVENTEENTH TERM
        ELSE
                SIGUV=0.0
                IPATH=17
        ENDIF

        RETURN 
        END

C***********************************************************************
C       THIS ROUTINE REDEFINES PARAMETERS INTO THE THE
C       THE TERMINOLOGY ASSOCIATED WITH DR. LII'S PAPER
C***********************************************************************
C       INPUT
C               R=NPO, RP=NTO, S=NPI, SP=NTI
C               P=R+S, PP=RP+SP, NETA=P+PP+1
C               A=1/PO, B=1/PI, AP=1/TO, BP=1/TI
C               MAXL----MAXIMUM LENGTH OF INVERSE POLYNOMIALS
C               ACURCY-THE ACURCY OF THE INVERSE POLYNOMIALS
C       REFER TO PUCOEF FOR THE REST OF THE DEFINITIONS
C***********************************************************************

      SUBROUTINE SHIFT(PI,NPI,PO,NPO,TI,NTI,TO,NTO,NP,MAXL,ACURCY,
     1       R,RP,S,SP,P,PP,NETA,A,NA,AP,NAP,B,NB,BP,NBP)

        INTEGER NPI,NPO,NTI,NTO,R,RP,S,SP,P,PP,NETA,NA,NAP,NB,NBP
        INTEGER NP,MAXL,I
        REAL PI(NP),PO(NP),TI(NP),TO(NP), ACURCY
        REAL A(0:MAXL), AP(0:MAXL), B(0:MAXL), BP(0:MAXL)  

        NA=0
        NB=0
        NAP=0
        NBP=0

C  CREATE THE INVERTED POLYNOMIALS AND SHIFT

        DO 10 I=0,MAXL
                A(I)=0.0
                B(I)=0.0
                AP(I)=0.0
                BP(I)=0.0
10      CONTINUE


        IF(NPO.GT.0) THEN
                CALL POLYIV(1,PO,NPO,NP,
     1                      A,MAXL,MAXL,
     2                      ACURCY,NA)
        ENDIF

        IF(NPI.GT.0) THEN
                CALL POLYIV(2,PI,NPI,NP,
     1                      B,MAXL,MAXL,
     2                      ACURCY,NB)
                DO 20 I=NPI+NB,NPI,-1
                       B(I)=B(I-NPI)
20              CONTINUE
                DO 30 I=0,NPI-1
                       B(I)=0
30              CONTINUE
                NB=NB+NPI
        ENDIF

        IF(NTO.GT.0) THEN
                CALL POLYIV(1,TO,NTO,NP,
     1                      AP,MAXL,MAXL,
     2                      ACURCY,NAP)
        ENDIF

        IF(NTI.GT.0) THEN
                CALL POLYIV(2,TI,NTI,NP,
     1                      BP,MAXL,MAXL,
     2                      ACURCY,NBP)
        DO 40 I=NTI+NBP,NTI,-1
               BP(I)=BP(I-NTI)
40      CONTINUE
        DO 50 I=0,NTI-1
               BP(I)=0
50      CONTINUE
        NBP=NBP+NTI
        ENDIF

C  PUT INTO FORM FOR USE WITH DR. LIIS ASYMPTOTIC COVARIANCE MATRIX
        R=NPO
        RP=NTO
        S=NPI
        SP=NTI
        P=R+S
        PP=RP+SP
        NETA=P+PP+1

        RETURN
        END

C00000000111111111122222222223333333333444444444455555555556666666666777
C23456789012345678901234567890123456789012345678901234567890123456789012
C00000000111111111122222222223333333333444444444455555555556666666666777
C23456789012345678901234567890123456789012345678901234567890123456789012
C
C  THESE ROUTINES ARE SIMPLE UTILITY ROUTINES USED THROUGHOUT THE
C  PROGRAM
C**********************************************************************
C        MATMUL--MULTIPLIES TWO MATRICES
C        TRNSPS--TRANSPOSES A MATRIX
C        XMEAN---MEAN OF A SEQUENCE
C       XSTDEV--STANDARD DEVIATION OF A SEQUENCE
C**********************************************************************
C
C        MATRIX MULTIPLICATION WITH C=A*B
C**********************************************************************
C
C        INPUT
C                A---------FIRST MATRIX
C                NRA,NCA---NUMBER OF ROWS AND COLUMNS OF A
C                B---------SECOND MATRIX
C                NRB,NCB---NUMBER OF ROWS AND COLOUNS OF B
C                NP--------DIMENSION OF THE MATRICES
C        OUTPUT
C                C---------PRODUCT OF THE ABOVE MATRICES
C
C        NOTE------IT IS ASSUMED THAT NCA=NRB
C                        NO CHECKING OF THIS IS DONE
C**********************************************************************

      SUBROUTINE MATMUL(A,NRA,NCA,B,NRB,NCB,C,NP)
        INTEGER NRA,NCA,NRB,NCB,NP,I,J,K
        REAL A(NP,NP), B(NP,NP), C(NP,NP)

        IF(NRB.EQ.NCA) THEN
        DO 30 I=1,NRA
                DO 20 J=1,NCB
                        C(I,J)=0.0
                        DO 10 K=1,NCA
                                C(I,J)=C(I,J)+A(I,K)*B(K,J)
10                        CONTINUE
20                CONTINUE
30      CONTINUE
        ENDIF

        RETURN
        END

C**********************************************************************
C          THIS SUBROUTINE CREATES THE TRANSPOSE OF A MATRIX
C**********************************************************************
C
C        INPUT
C                A---------FIRST MATRIX
C                NR,NC---NUMBER OF ROWS AND COLUMNS OF A
C                NP--------DIMENSION OF THE MATRICES
C        OUTPUT
C                B---------SECOND MATRIX
C**********************************************************************

      SUBROUTINE TRNSPS(A,NR,NC,B,NP)
        INTEGER NR,NC,NP,I,J
        REAL A(NP,NP),B(NP,NP)

        DO 20 I=1,NR
                DO 10,J=1,NC
                        B(J,I)=A(I,J)
10                CONTINUE
20        CONTINUE


        RETURN
        END
C**********************************************************************
C          THIS FUNCTION FINDS THE MEAN OF A TRUNCATED SEQUENCE
C**********************************************************************
C
C        INPUT
C                X---------SEQUENCE VECTOR
C                NX--------SEQUENCE DIMENSION
C                                OR SIMPLY UP TO THE NUMBER OF VALUES
C                                IN THE SEQUENCE
C                IQ--------FRONT AND BACK TRUNCATION LENGTH
C**********************************************************************

      FUNCTION XMEAN(X,NX,IQ)
        INTEGER NX,IQ,I
        REAL X(NX),SUM,TOT,XMEAN

        SUM=0
        TOT=NX-2*IQ
        DO 10 I=IQ+1,NX-IQ
                SUM=SUM+X(I)
10        CONTINUE
        XMEAN=SUM/TOT

        RETURN
        END

C**********************************************************************
C          THIS FUNCTION FINDS THE STANDARD DEVIATION OF A 
C        TRUNCATED SEQUENCE
C**********************************************************************
C
C        INPUT
C                X---------SEQUENCE VECTOR
C                NX--------SEQUENCE DIMENSION
C                                OR SIMPLY UP TO THE NUMBER OF VALUES
C                                IN THE SEQUENCE
C                IQ--------FRONT AND BACK TRUNCATION LENGTH
C**********************************************************************


      FUNCTION XSTDEV(X,NX,IQ)
        INTEGER NX,IQ,I
        REAL X(NX),SUM2,TOT,XM,XMEAN,XVAR,XSTDEV

        SUM2=0.0
        TOT=NX-2*IQ
        DO 10 I=IQ+1,NX-IQ
                SUM2=SUM2+X(I)*X(I)
10        CONTINUE

        XM=XMEAN(X,NX,IQ)
        XVAR=(SUM2-TOT*XM*XM)/(TOT-1)
        XSTDEV=SQRT(XVAR)
        RETURN        
        END

C00000000111111111122222222223333333333444444444455555555556666666666777
C23456789012345678901234567890123456789012345678901234567890123456789012
