PROGRAM BackProp
! Copyright 2007 Matthew Norton
!
!    This program is free software; you can redistribute it and/or modify
!    it under the terms of the GNU General Public License as published by
!    the Free Software Foundation; either version 2 of the License, or
!    (at your option) any later version.
!
!    This program is distributed in the hope that it will be useful,
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!    GNU General Public License for more details.
!
!    You should have received a copy of the GNU General Public License
!    along with this program; if not, write to the Free Software
!    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
!
! NOTE: This is in FORTRAN 90 and you will need a FORTRAN 90 complier to run this code
!---------------------------------------------------------------------------
! This program uses backpropagation to train an artificial neural network
! to learn the Dow Jones Industrial Average for the first 200 months of 
! data. These weights can be used to attempt to predict future values of the 
! market
!
! Programmed by:
!     Matthew Norton
!
!----------------------------------------------------------------------------
IMPLICIT NONE
INTEGER :: I, J, M = 0, N = 0, Count = 0, X, K, ECount = 1, B, MaxCount = 1000000, DWCount, DVWCount
REAL*8 :: Alpha, G, FPrimeGk
REAL*8, DIMENSION(:), ALLOCATABLE :: DJIA
REAL*8, DIMENSION(10) :: Zj = 0.0, Hj = 0.0
REAL*8 :: EpsilonK = 0.0, Gk = 0.0
REAL*8, DIMENSION(0:10,6) :: Weights, DWeights
REAL*8, DIMENSION(0:6) :: VWeights, DVWeights, DeltaJ
REAL*8, DIMENSION(6) :: FPRimeZj
REAL*8, DIMENSION(200,10) :: Train = 0.0
REAL*8, DIMENSION(200) :: Tar, Yk, Error, RError
 CALL SetupDJIA()
 CALL Setup()
      DO J = 1, 200
           DO X = 1, 10 
                N = X + M ! Increment counter for DJIA
                Train(J,X) = DJIA(N) ! Train(J,I) = DJIA(I+M)
!                PRINT *, J, X, N 
           END DO
      Tar(J) = DJIA(N + 1) ! Target output = DJIA(I+M+1)
      !PRINT *, Tar(J)
      M = M + 1 ! Increment M
      END DO
 Epoch: DO
      WeightUpdate: DO B = 1, 200
               DO J = 1, 6
                   Zj(J) = Weights(0,J) + (Weights(1,J) * Train(B,1)) +&
                   &(Weights(2,J) * Train(B,2)) + (Weights(3,J) * Train(B,3)) +&
                   &(Weights(4,J) * Train(B,4)) + (Weights(5,J) * Train(B,5)) +&
                   &(Weights(6,J) * Train(B,6)) + (Weights(7,J) * Train(B,7)) +&
                   &(Weights(8,J) * Train(B,8)) + (Weights(9,J) * Train(B,9)) +&
                   &(Weights(10,J) * Train(B,10)) 
                   Hj(J) = (( 1.0 + DTANH(Zj(J))) / 2.0)
               END DO
          Gk = VWeights(0)
          DO J = 1, 6
              Gk = Gk + (VWeights(J) * Hj(J))
          END DO
          Yk(B) = (( 1.0 + DTANH(Gk)) / 2.0)
          FprimeGk = 1.0 / (2 * DCOSH(Gk) * DCOSH(Gk))
          EpsilonK = (Tar(B) - Yk(B)) * FPrimeGk
          DVWeights(0) = Alpha * EpsilonK
          DO J = 1, 6
               DVWeights(J) = Alpha * EpsilonK * Hj(J)
          END DO
          DO J = 1, 6
               FPrimeZj(J) = 1.0 / (2 * DCOSH(Zj(J)) * DCOSH(Zj(J)))
               DeltaJ(J) = EpsilonK * VWeights(J) * FPrimeZj(J)
          END DO
          DO I = 1, 10
               DO J = 1, 6
                    DWeights(I,J) = Alpha * DeltaJ(J) * Tar(B)
               END DO
          END DO
          DO J = 1, 6
               DWeights(0,J) = Alpha * DeltaJ(J)
          END DO
          DO I = 0, 10
               DO J = 1, 6
                    Weights(I,J) = Weights(I,J) + DWeights(I,J)
               END DO
          END DO
          DO J = 0, 6
               VWeights(J) = VWeights(J) + DVWeights(J)
          END DO
     END DO WeightUpdate
     DO B = 1, 200
          DO J = 1, 6
               Zj(J) = Weights(0,J) + (Weights(1,J) * Train(B,1)) +&
               &(Weights(2,J) * Train(B,2)) + (Weights(3,J) * Train(B,3)) +&
               &(Weights(4,J) * Train(B,4)) + (Weights(5,J) * Train(B,5)) +&
               &(Weights(6,J) * Train(B,6)) + (Weights(7,J) * Train(B,7)) +&
               &(Weights(8,J) * Train(B,8)) + (Weights(9,J) * Train(B,9)) +&
               &(Weights(10,J) * Train(B,10)) 
               Hj(J) = (( 1.0 + DTANH(Zj(J))) / 2.0)
          END DO
          Gk = VWeights(0)
          DO J = 1, 6
               Gk = Gk + (VWeights(J) * Hj(J))
          END DO
          Yk(B) = (( 1.0 + DTANH(Gk)) / 2.0)
     END DO
     DO B = 1, 200
          Rerror(B) = ABS(Yk(B) - Tar(B)) / Tar(B)
     END DO
     PRINT *, ECount
     ECount = Ecount + 1
     IF(Ecount > MaxCount) THEN
          PRINT *, "Failed to converge after", ECount, "Epochs"
          EXIT Epoch
     END IF
     DO B = 1, 200
          IF (RError(B) > 0.05) THEN
               CALL Reset
               CYCLE Epoch
          END IF
     END DO
     DO B = 1, 200
          Error(B) = 0.5 * ((Yk(B) - Tar(B)) * (Yk(B) - Tar(B)))
          Rerror(B) = ABS(Yk(B) - Tar(B)) / Tar(B)
          WRITE(20, 200) B, Yk(B), Tar(B), Error(B), RError(B)
     END DO

     DO I = 0, 10
          DO J = 1, 6
               WRITE(30, 300, ADVANCE = "NO") Weights(I,J)
          END DO
     END DO
     DO J = 0, 6
          WRITE(40, 300, ADVANCE = "NO") VWeights(J)
     END DO
     PRINT *, "The computer has learned the stock market"
     EXIT Epoch
 END DO Epoch     

200 FORMAT(1X, I5, 3(1X, F15.13), F15.10)
300 FORMAT(1X, F15.12)

! ------------------------------------------------------------------------------
 CONTAINS
! ------------------------------------------------------------------------------

SUBROUTINE SetupDJIA
INTEGER :: I = 1, ErrorStat
REAL :: Dummy ! Dummy variable used for reading in array
OPEN(UNIT = 10, FILE = "DJI.txt", STATUS = "OLD", IOSTAT = ErrorStat) !input file
IF (ErrorStat /= 0) STOP "*** ERROR IN OPENING DJI.txt ***" ! If error, stop
FindSize: DO ! Find file length 
     READ (10,200) Dummy ! read in line
     READ (10,200, IOSTAT = ErrorStat, ADVANCE = "NO") Dummy ! see if next line is there
     IF (ErrorStat /= 0) THEN ! If end of file, exit loop
          EXIT FindSize
     END IF
     Count = Count + 1 ! Increment count
END DO FindSize
REWIND(10) ! go back to begining of file
ALLOCATE(DJIA(Count)) ! Allocate space for DJIA array
DO I = 1, Count ! Do from I to end of file
     READ(10,200) DJIA(I) ! read in values for DJIA from file
END DO
200 FORMAT(T2, F8.7)
END SUBROUTINE

SUBROUTINE Setup
REAL*8 :: R1, R2, Sigma = 0.0, Sum = 0.0
REAL*8, DIMENSION(66) :: Sums = 0.0
REAL, DIMENSION(6) :: VSums = 0.0
INTEGER :: Sign = 1, I
OPEN(UNIT = 20, File = "BackOut.txt", Status = "UNKNOWN")
OPEN(UNIT = 30, File = "Weights.txt", Status = "UNKNOWN")
OPEN(UNIT = 40, FILE = "VWeights.txt", STATUS = "UNKNOWN")
OPEN(UNIT = 50, FILE = "Dweights.txt", STATUS = "UNKNOWN")
OPEN(UNIT = 60, FILE = "DVWeights.txt", STATUS = "UNKNOWN")
 CALL Seed()
 Sigma = 0.7 * 6.0 ** (1.0 / 10.0)
 Alpha = 0.05
DO I = 1, 10
     DO J = 1, 6
     Sign = 1
     CALL RANDOM_NUMBER(R1)
     CALL RANDOM_NUMBER(R2)
     IF (R2 < 0.5) Sign = -1
     DO
          IF (R1 <= (1.0 / SQRT(11.0))) THEN
               EXIT
          END IF
          CALL RANDOM_NUMBER(R1)
     END DO
     Weights(I,J) = Sign * R1
     END DO
END DO
DO I = 1, 10
     DO J = 1, 6
          Sum = Sum + Weights(I,J) * Weights(I,J)
     END DO
     Sums(I) = SQRT(Sum)
     Sum = 0.0
END DO
DO I = 1, 10
     DO J = 1, 6
          Weights(I,J) = (Sigma * Weights(I,J)) / Sums(I)
     END DO
END DO
DO J = 1, 6
     Sign = 1
     CALL RANDOM_NUMBER(R1)
     CALL RANDOM_NUMBER(R2)
     IF (R2 < 0.5) Sign = -1
     DO
          IF (R1 <= Sigma) THEN
               EXIT
          END IF
          CALL RANDOM_NUMBER(R1)
     END DO
     Weights(0,J) = Sign * R1
     END DO
Sum = 0.0
DO I = 1, 6
     Sign = 1
     CALL RANDOM_NUMBER(R1)
     CALL RANDOM_NUMBER(R2)
     IF (R2 < 0.5) Sign = -1
     DO
          IF (R1 <= (1.0 / SQRT(7.0))) THEN
               EXIT
          END IF
          CALL RANDOM_NUMBER(R1)
     END DO
     VWeights(I) = Sign * R1
END DO
DO I = 1, 6
     Sum = Sum + VWeights(I) * VWeights(I)
     VSums(I) = SQRT(Sum)
     Sum = 0.0
END DO
DO I = 1, 6
     VWeights(I) = (Sigma * VWeights(I)) / Sums(I)
END DO
 Sign = 1
 CALL RANDOM_NUMBER(R1)
 CALL RANDOM_NUMBER(R2)
 IF (R2 < 0.5) Sign = -1
 DO
     IF (R1 <= Sigma) THEN
          EXIT
     END IF
     CALL RANDOM_NUMBER(R1)
 END DO
 VWeights(0) = Sign * R1
END SUBROUTINE

! ------------------------------------------------------------------------------

SUBROUTINE Seed
INTEGER :: Now(8)
 CALL date_and_time(values=Now)
 CALL random_seed(put=Now(5:8))
END SUBROUTINE

! ------------------------------------------------------------------------------

SUBROUTINE Reset
DWeights = 0.0
DVWeights = 0.0
Error = 0.0
Zj = 0.0
Hj = 0.0
EpsilonK = 0.0
Gk = 0.0
Yk = 0.0
END SUBROUTINE
END PROGRAM