PROGRAM GENDAT ! ------------------------------------------------
      PARAMETER (N=100,M=5,K=6,NCUE=1)
C     N=NO. OF OBSERVATIONS, M=NO. OF EQUATIONS (= NO. OF ENDOGENOUS
C     VARIABLES, Y) AND K=NO. OF EXOGENOUS VARIABLES (X) INCLUDING A
C     UNITARY VARIABLE ASSOCIATED WITH THE INTERCEPT IN THE KTH COLUMN.
C     IF NCUE=1 THEN U WITH N(0,SDU) IS DIRECTLY GENERATED & ADDED TO Y
C     IF NCUE=0, U IS -E*INV(A), E WITH N(0,SDE). THEN U IS ADDED TO Y
      IMPLICIT DOUBLE PRECISION  (A-H,O-Z)
      PARAMETER (XLOW=0,XRANGE=20, SDU=0.0010D0,SDE=1.818D-07)
C     XLOW=LOWEST LIMIT ON X(I,J); XRANGE=RANGE OF X(I,J) FROM XLOW TO
C     XHIGH, THE UPPER LIMIT OF X(I,J). SDU=STANDARD DEVIATION OF NORMAL
C     ERROR U WITH N(0, SDU), SDE OF E WITH N(0,SDE).
      CHARACTER *30 OFIL
      PARAMETER (OFIL='TSLDAT.TXT')
      COMMON /RNDM/IU,IV
      DIMENSION Y(N,M),X(N,K),A(M,M),B(K,M),P(K,M),E(N,M),U(N,M)
C     --------------   Structural Parameters   -----------------------
      DATA ((A(I,J),I=1,5),J=1,5)/-1,7,0,-6,0,  3,-1,5,0,0,  0,0,-1,3,0,
     & 6,0,0,-1,-3,  -11,0,9,0,-1/
      DATA ((B(I,J),I=1,6),J=1,5)/0,5,0,-7,0,60,  3,0,-5,0,0,20,
     & 0,2,0,0,0,9,  0,4,0,0,-3,-8,  0,0,0,6,0, -11/
C     -----------------------------------------------------------------
      WRITE(*,*)'FEED A 4-DIGIT ODD NON-ZERO RANDOM NUMBER SEED'
      READ(*,*) IU
C     PRINT TRANSPOSE(A) AND TRANSPOSE(B)
C     -----------------------------------------------------------------
      WRITE(*,*)'TRANSPOSED A MATRIX'
      DO J=1,M
      WRITE(*,10)(A(I,J),I=1,M)
      ENDDO
      WRITE(*,*)' '
      WRITE(*,*)'TRANSPOSED B MATRIX'
      DO J=1,M
      WRITE(*,10)(B(I,J),I=1,K)
      ENDDO
      WRITE(*,*)' '
   10 FORMAT(6F7.0)
C     -----------------------------------------------------------------
C     GENERATE X(N,K)
      DO I=1,N
      DO J=1,K-1
      CALL RANDOM(RAND)
      X(I,J)=XLOW+RAND*XRANGE
      ENDDO
      ENDDO
      DO I=1,N
      X(I,K)=1.D0
      ENDDO
C     INVERT A MATRIX AND FIND P=-B*INV(A)
      CALL MINV(A,M,D)
      WRITE(*,*)'DETERMINANT=',D
      DO I=1,K
      DO J=1,M
      P(I,J)=0.D0
      DO JJ=1,M
      P(I,J)=P(I,J)+B(I,JJ)*A(JJ,J)
      ENDDO
      ENDDO
      ENDDO
C     PRINT P MATRIX
      WRITE(*,*)'REDUCED FORM COEEFIENTS OR  P MATRIX ---------------'
      DO I=1,K
      WRITE(*,1)(P(I,J),J=1,M)
      ENDDO
    1 FORMAT(6F12.4)
C     OBTAIN Y= XP
      DO I=1,N
      DO J=1,M
      Y(I,J)=0.D0
      DO JJ=1,K
      Y(I,J)=Y(I,J)+X(I,JJ)*P(JJ,J)
      ENDDO
      ENDDO
      ENDDO
C     PRINT Y (WITHOUT ERROR)
      WRITE(*,*)'-------Y WITHOUT ERROR ------------------------------'
      DO I=1,N
      WRITE(*,1)(Y(I,J),J=1,M)
      ENDDO
C     --------------- ADDING U WITH N(0,SDU) TO Y DIRECTLY -----------
      IF(NCUE.EQ.1) THEN
C     ADD NORMAL ERRORS TO Y
      DO J=1,M
      DO I=1,N
      CALL NORMAL(RAND) ! GENERATE NORMALLY DISTRIBUTED ERRORS N(0,SDE)
      Y(I,J)=Y(I,J)+RAND*SDU ! ADD ERROR TO Y
      ENDDO
      ENDDO
      ENDIF
C     --- GENERATING U =  -E INV(A) AND ADDING TO Y ------------------
      IF(NCUE.EQ.0) THEN
      DO I=1,N
      DO J=1,M
      CALL NORMAL(RAND)
      E(I,J)=RAND*SDE
      ENDDO
      ENDDO
      DO I=1,N
      DO J=1,M
      U(I,J)=0.D0
      DO JJ=1,M
      U(I,J)=U(I,J)-E(I,JJ)*A(JJ,J)
      ENDDO
      ENDDO
      ENDDO
      DO I=1,N
      DO J=1,M
      Y(I,J)=Y(I,J)+U(I,J)
      ENDDO
      ENDDO
      ENDIF
C     -----------------------------------------------------------------
C     PRINT Y (WITH ERROR ADDED)
      WRITE(*,*)'-------Y WITH ERROR ------------------------------'
      DO I=1,N
      WRITE(*,1)(Y(I,J),J=1,M)
      ENDDO
      OPEN(7,FILE=OFIL)
      DO I=1,N
      WRITE(7,1)(Y(I,J),J=1,M)
      ENDDO
      DO I=1,N
      WRITE(7,1)(X(I,J),J=1,K)
      ENDDO
      WRITE(7,*)'P MATRIX -------------------------------------------'
      DO J=1,K
      WRITE(7,1)(P(J,JJ),JJ=1,M)
      ENDDO
      CLOSE(7)
      END
C     -----------------------------------------------------------------
C     RANDOM NUMBER GENERATOR (UNIFORM BETWEEN 0 AND 1 - BOTH EXCLUSIVE)
      SUBROUTINE RANDOM(RAND)
       IMPLICIT DOUBLE PRECISION  (A-H,O-Z)
       COMMON /RNDM/IU,IV
       IV=IU*65539
       IF(IV.LT.0) THEN
       IV=IV+2147483647+1
       ENDIF
       RAND=IV
       IU=IV
       RAND=RAND*0.4656613E-09
       RETURN
       END
C     -----------------------------------------------------------------
      SUBROUTINE NORMAL(R)
C     PROGRAM TO GENERATE N(0,1) FROM RECTANGULAR RANDOM NUMBERS
C     IT USES VARIATE TRANSFORMATION FOR THIS PURPOSE.
C     -----------------------------------------------------------------
C     IF U1 AND U2 ARE UNIFORMLY DISTRIBUTED RANDOM NUMBERS (0,1),
C     THEN X=[(-2*LN(U1))**.5]*(COS(2*PI*U2) IS N(0,1)
C     PI = 4*ARCTAN(1.0)= 3.1415926535897932384626433832795
C     2*PI = 6.283185307179586476925286766559
C     -----------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /RNDM/IU,IV
C     -----------------------------------------------------------------
      CALL RANDOM(RAND) ! INVOKES RANDOM TO GENERATE UNIFORM RAND [0, 1]
      U1=RAND ! U1 IS UNIFORMLY DISTRIBUTED [0, 1]
      CALL RANDOM(RAND) ! INVOKES RANDOM TO GENERATE UNIFORM RAND [0, 1]
      U2=RAND ! U1 IS UNIFORMLY DISTRIBUTED [0, 1]
      R=DSQRT(-2.D0*DLOG(U1))
      R=R*DCOS(U2*6.283185307179586476925286766559D00)
C     R=R*DCOS(U2*6.28318530718D00)
      RETURN
      END
C     -----------------------------------------------------------------
C     SUBROUTINE FOR MATRIX INVERSION
      SUBROUTINE MINV(A,N,D)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(N,N)
      U=1.D0
      D=U
      DO I=1,N
      D=D*A(I,I)
      A(I,I)=U/A(I,I)
         DO J=1,N
         IF(I.NE.J) A(J,I)=A(J,I)*A(I,I)
         ENDDO
         DO J=1,N
         DO K=1,N
         IF(I.NE.J.AND.K.NE.I) A(J,K)=A(J,K)-A(J,I)*A(I,K)
         ENDDO
         ENDDO
         DO J=1,N
         IF(J.NE.I) A(I,J)= -A(I,J)*A(I,I)
         ENDDO
      ENDDO
      RETURN
      END