PROGRAM VARCOMPUTE
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      INTEGER NIT(0:249,0:4), NKT(0:499,0:4)
      INTEGER NIK(2000),IIK(2000),KIK(2000)
      INTEGER T, LEVELS(6)
      CHARACTER*15 NAME(6)
      DIMENSION A(5,6),TADJ(5),XXX(11),xl(6)
      real*4 roa, r(8192)
      COMMON /KVCOM / KV(8192,4)
      common /rcom/ r
      COMMON /XCOM/ XMA,XMB,XMC,XMD,XME,XAB,XAC,XBA,XBC,XCA,XCB
      COMMON /LEVCOM/ LEVA,LEVB,LEVC,LEVD,LEVE,LEVERROR
      EQUIVALENCE (XMA,XXX(1))
      EQUIVALENCE (LEVA,LEVELS(1))
C
      DATA NIT/1250*0/, NKT/2500*0/, NIK/2000*0/
      DATA IIK/2000*0/, KIK/2000*0/
      
C
      DATA SIGE2 /102.66/
      DATA TADJ /347423.96,340475.23,548.91,537841.56,
     1  1294372.15/
      NAME(1) = ‘INDUSTRY’
      NAME(2) = ‘CORPORATION’
      NAME(3) = ‘TIME’
      NAME(4) = ‘TIME*INDUSTRY’
      NAME(5) = ‘SBU’
      NAME(6)=’ERROR’
      LEVERROR = 0

      DO 50 I=1,11
50    XXX(I) = 0.0

C     TA HOLDS THE CALCULATED VALUES OF T – TM
C     TADJ(1) = TA – TM
C     TADJ(2) = TB – TM
C     TADJ(3) = TC – TM
C     TADJ(4) = TD – TM
C     TADJ(5) = TE – TM
C
C
      OPEN(9,FILE=’VARCOMPUTE.LOG’)
      CLOSE (9,STATUS=’DELETE’)
      OPEN(9,FILE=’VARCOMPUTE.LOG’)
      WRITE(9,10)
10    FORMAT(‘ VARIANCE COMPONENTS RESULTS’//)
        print ‘(” Opening lbd file”)’
      OPEN(8,FILE=’d:\data\LBD1.dat’,STATUS=’OLD’)
      NOB = 0
C
C     READ INTO KV(.,1) IYEAR (C)
C     READ INTO KV(.,2) NCORP (B)
C     READ INTO KV(.,3) IND4  (A)
C     READ INTO KV(.,4) NSBU  (E)
C
      tz=0.0
      tm=0.0
      ta=0.0
      tb=0.0
      tc=0.0
      td=0.0
      tf=0.0
200   READ(8,202,END=220) IYEAR, NCORP, IND4, roa
202   FORMAT(I2,I4,I5,f8.3)
      NOB=NOB+1
      tz = tz+roa*roa
      tm = tm+roa
      r(nob)=roa
      KV(NOB,1)= IYEAR
      KV(NOB,2)= NCORP
      KV(NOB,3)= IND4
      KV(NOB,4)= 10000*NCORP + IND4
            GOTO 200
220   CLOSE (8)
      PRINT ‘(I5,” OBSERVATIONS”)’,NOB
      tm = tm/dble(nob)
C
      CALL LEVELER(NOB,LEVA,KV(1,3))
      CALL LEVELER(NOB,LEVB,KV(1,2))
      CALL LEVELER(NOB,LEVC,KV(1,1))
      CALL LEVELER(NOB,LEVE,KV(1,4))
      LEVD = 0
      DO 250 N=1,NOB
        T = KV(N,1)
        I = KV(N,3)
        K = KV(N,2)
        J = KV(N,4)
        NIT(I,T) = NIT(I,T) + 1
        IF(NIT(I,T).EQ.1) LEVD = LEVD + 1
        NIT(I,0) = NIT(I,0) + 1
        NIT(0,T) = NIT(0,T) + 1
        NKT(K,T) = NKT(K,T) + 1
        NKT(0,T) = NKT(0,T) + 1
        NKT(K,0) = NKT(K,0) + 1
        IF(NIK(J).EQ.0) THEN
          NIK(J) = 1
          IIK(J) = I
          KIK(J) = K
        ELSE
          NIK(J) = NIK(J) + 1
          IF(IIK(J).NE.I.OR.KIK(J).NE.K) THEN
            PRINT ‘(” DATA CONFUSION PROBLEM”)’
            GOTO 1000
          ENDIF
        ENDIF
250   CONTINUE
C       nit holds obs in (i,t)
c       nkt holds obs in (k,t)
c       nik holds obs in bus unit (j)
c           iik tells i associated with j
c           kik tells k associated with j
c    
 
      DO 275 J = 1,LEVE
        FNIKD2 = NIK(J)*NIK(J)
        XME = XME + FNIKD2/NOB
        CALL CONDADD(XAB,FNIKD2,NIT(IIK(J),0))
        CALL CONDADD(XBA,FNIKD2,NKT(KIK(J),0))
275   CONTINUE
C
      DO 300 I = 1,LEVA
        FNIDD = NIT(I,0)
        XMA = XMA + FNIDD*FNIDD/NOB
        DO 300 T = 1,LEVC
          FNIDT2 = NIT(I,T)*NIT(I,T)
        XMD = XMD + FNIDT2/NOB
          CALL CONDADD(XAC,FNIDT2,NIT(I,0))
          CALL CONDADD(XCA,FNIDT2,NIT(0,T))
300   CONTINUE
C
      DO 350 K = 1,LEVB
        FNDKD = NKT(K,0)
        XMB = XMB + FNDKD*FNDKD/NOB
        DO 350 T=1,LEVC
          FNDKT2 = NKT(K,T)*NKT(K,T)
          CALL CONDADD(XBC,FNDKT2,NKT(K,0))
          CALL CONDADD(XCB,FNDKT2,NKT(0,T))
350   CONTINUE
      DO 400 T=1,LEVC
        FNDDT2 = NKT(0,T)*NKT(0,T)
400   XMC = XMC + FNDDT2/NOB
C
      A(1,1) = NOB – XMA
      A(1,2) = XAB – XMB
      A(1,3) = XAC – XMC
      A(1,4) = XAC – XMD
      A(1,5) = XAB – XME
      A(1,6) = TADJ(1) – (LEVA-1)*SIGE2
      A(2,1) = XBA – XMA
      A(2,2) = NOB – XMB
      A(2,3) = XBC – XMC
      A(2,4) = LEVB – XMD
      A(2,5) = XBA – XME
      A(2,6) = TADJ(2) – (LEVB-1)*SIGE2
      A(3,1) = XCA – XMA
      A(3,2) = XCB – XMB
      A(3,3) = NOB – XMC
      A(3,4) = XCA – XMD
      A(3,5) = LEVC – XME
      A(3,6) = TADJ(3) – (LEVC-1)*SIGE2
      A(4,1) = NOB – XMA
      A(4,2) = LEVD – XMB
      A(4,3) = NOB – XMC
      A(4,4) = NOB – XMD
      A(4,5) = LEVD – XME
      A(4,6) = TADJ(4) – (LEVD-1)*SIGE2
      A(5,1) = NOB – XMA
      A(5,2) = NOB – XMB
      A(5,3) = LEVE – XMC
      A(5,4) = LEVE – XMD
      A(5,5) = NOB – XME
      A(5,6) = TADJ(5) – (LEVE-1)*SIGE2
501   format(7f20.9)
        xl(1) = leva
        xl(2) = levb
        xl(3) = levc
        xl(4) = levd
        xl(5) = leve
        write(9,501) nob-xma,nob-xmb,nob-xmc,nob-xmd,nob-xme,nob-1.0,0.0
        do 510 i = 1,5
      write(9,501) (a(i,j),j=1,6),xl(i)-1,tadj(i)
510     continue
      CALL SOLVER(5,6,A)
601   FORMAT(‘SOURCE’,T17,’LEVELS’,T34,’ADJ SSQ R’,T50,
     1 ‘VARIANCE COMPONENT’)
602   FORMAT(A15,I7,2F20.6)
603   FORMAT(‘ERROR’,T43,F20.6)
      WRITE(9,601)
      DO 600 I=1,5
600   WRITE(9,602) NAME(I),LEVELS(I),TADJ(I),A(I,6)
      WRITE(9,603) SIGE2
      WRITE(*,601)
      DO 610 I=1,5
610   WRITE(*,602) NAME(I),LEVELS(I),TADJ(I),A(I,6)
      WRITE(*,603) SIGE2
1000  CLOSE(8)
      CLOSE (9)
      END
C
      SUBROUTINE CONDADD(X,F,N)
      DOUBLE PRECISION X,F
      IF(N.GT.0) X = X + F/N
      RETURN
      END
C     SUBROUTINE LEVELER
C
C     THIS SUBROUTINE ACCEPTS AN INTEGER VECTOR
C     CODES(.) AND REPLACES EACH UNIQUE CODE WITH
C     ITS ORDER NUMBER (LEVEL) ON AN ASCENDING SORT.
C     THE NUMBER OF LEVELS (LARGEST ORDER NUMBER) IS
C     RETURNED IN NLEV AND THE VECTOR IW(.) IS USED AS
C     WORKSAPCE.  CODES AND IW ARE ADJ TO SIZE N.
C
      SUBROUTINE LEVELER(N,NLEV,CODES)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      INTEGER CODES(N)
      INTEGER IW(8192)
      CALL INDEXX(N,CODES,IW)
      NLEV = 1
      LAST = CODES(IW(1))
      DO 100 I =1,N
        IF(LAST.NE.CODES(IW(I))) THEN
          NLEV = NLEV + 1
          LAST = CODES(IW(I))
        ENDIF
        CODES(IW(I)) = NLEV
100   CONTINUE
      RETURN
      END
C
C     SUBROUTINE INDEXX
C     THIS SUBROUTINE SORTS THE INTEGER VECTOR A(.)
C     VIA AN INDEX VECTOR INDX(.).  THE ADJUSTED SIZES
C     OF THESE VECTORS IS N.
C     A HEAPSORT METHOD IS USED.  THE CONTENTS OF A(.)
C     ARE NOT ALTERED.
C
      SUBROUTINE INDEXX(N,A,INDX)
      INTEGER A(N), INDX(N), Q
      DO 100 J =1,N
100   INDX(J) = J
      L = N/2 + 1
      IR = N
110   CONTINUE
        IF(L.GT.1) THEN
          L = L-1
          INDXT = INDX(L)
          Q = A(INDXT)
        ELSE
          INDXT = INDX(IR)
          Q = A(INDXT)
          INDX(IR) = INDX(1)
          IR = IR-1
          IF(IR.EQ.1) THEN
            INDX(1) = INDXT
            RETURN
          ENDIF
        ENDIF
        I=L
        J=L+L
120     IF(J.LE.IR) THEN
          IF(J.LT.IR) THEN
            IF(A(INDX(J)).LT.A(INDX(J+1))) J=J+1
          ENDIF
          IF(Q.LT.A(INDX(J))) THEN
            INDX(I)=INDX(J)
            I=J
            J=J+J
          ELSE
            J=IR+1
          ENDIF
        GOTO 120
        ENDIF
        INDX(I) = INDXT
      GOTO 110
      END
C
      SUBROUTINE SOLVER(N,M,A)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION A(N,M)
C
C     DO EACH COLUMN IN TURN
C
      DO 500 JP = 1,N
C
C       SEARCH FOR BIGGEST ELEMENT TO PIVOT ON
        IP = 0
        BIG = 0.0
        DO 100 I = JP,N
          IF(DABS(A(I,JP)).GT.BIG) THEN
            BIG = dabs(A(I,JP))
            IP = I
          ENDIF
100     CONTINUE
        IF(IP.EQ.0) THEN
          STOP
        ENDIF
C
C       IF NOT ON DIAGONAL, EXCHANGE ROWS
120     IF(IP.NE.JP) THEN
          DO 140 J=1,M
            X = A(IP,J)
            A(IP,J) = A(JP,J)
            A(JP,J) = X
140       CONTINUE
          IP = JP
        ENDIF
C
C       MAKE THE PIVOT UNITY
C
      Q = A(IP,JP)
        DO 200 J=JP,M
200     A(IP,J) = A(IP,J)/Q
C
C       CLEAR REST OF COLUMN
C
        DO 250 I = 1, N
          IF(I.EQ.IP) GOTO 250
          Q = A(I,JP)
          A(I,JP) = 0.0
          DO 220 J=JP+1,M
220       A(I,J) = A(I,J) – Q*A(IP,J)
250     CONTINUE
500   CONTINUE
      RETURN
      END