program regression
implicit none

real                      :: a
real                      :: b1
real                      :: b2
real                      :: fTest
integer                   :: i
integer                   :: k
integer                   :: n
real                      :: rSquared
real                      :: ssReg
real                      :: ssRes
real                      :: sumX1
real                      :: sumX1X1
real                      :: sumX1X1_norm
real                      :: sumX1X2
real                      :: sumX1X2_norm
real                      :: sumX1Y
real                      :: sumX1Y_norm
real                      :: sumX2
real                      :: sumX2X2
real                      :: sumX2X2_norm
real                      :: sumX2Y
real                      :: sumX2Y_norm
real                      :: sumY
real                      :: sumYY
real                      :: sumYY_norm
real,         allocatable :: x1(:)
real,         allocatable :: x2(:)
real,         allocatable :: y(:)
character(3), allocatable :: yName(:)


! Parameters
k = 2
n = 30


! Allocations
allocate(yName(n), y(n), x1(n), x2(n))

sumX1   = 0
sumX1X1 = 0
sumX1X2 = 0
sumX1Y  = 0
sumX2   = 0
sumX2X2 = 0
sumX2Y  = 0
sumY    = 0
sumYY   = 0


! Data
yName = [ "WSH", "PIT", "CHI", "CBJ", "MIN", "ANA", "MTL", "EDM", "NYR", "STL", &
        & "SJS", "OTT", "BOS", "TOR", "CGY", "NSH", "NYI", "TBL", "PHI", "WPG", &
        & "CAR", "LAK", "FLA", "DAL", "DET", "BUF", "NJD", "ARI", "VAN", "COL"  ]

y     = [ 55, 50, 50, 50, 49, 46, 47, 47, 48, 46, 46, 44, 44, 40, 45, &
        & 41, 41, 42, 39, 40, 36, 39, 35, 34, 33, 33, 28, 30, 30, 22  ]

x1    = [ 263, 282, 244, 249, 266, 223, 226, 247, 256, 235, 221, 212, 234, 251, 226, &
        & 240, 241, 234, 219, 249, 215, 201, 210, 223, 207, 201, 183, 197, 182, 166  ]

x2    = [ 182, 234, 213, 195, 208, 200, 200, 212, 220, 218, 201, 214, 212, 242, 221, &
        & 224, 242, 227, 236, 256, 236, 205, 237, 262, 244, 237, 244, 260, 243, 278  ]


! Preliminary calculations
do i = 1, n
   sumY    = sumY    + y(i)
   sumX1   = sumX1   + x1(i)
   sumX2   = sumX2   + x2(i)
   sumYY   = sumYY   + y(i)  * y(i)
   sumX1X1 = sumX1X1 + x1(i) * x1(i)
   sumX2X2 = sumX2X2 + x2(i) * x2(i)
   sumX1Y  = sumX1Y  + x1(i) * y(i)
   sumX2Y  = sumX2Y  + x2(i) * y(i)
   sumX1X2 = sumX1X2 + x1(i) * x2(i)
enddo

sumYY_norm   = sumYY   - (sumY)**2       / n
sumX1X1_norm = sumX1X1 - (sumX1)**2      / n
sumX2X2_norm = sumX2X2 - (sumX2)**2      / n
sumX1Y_norm  = sumX1Y  - (sumX1 * sumY)  / n
sumX2Y_norm  = sumX2Y  - (sumX2 * sumY)  / n
sumX1X2_norm = sumX1X2 - (sumX1 * sumX2) / n


! Calculate regression slopes and intercept
b1 = (sumX2X2_norm * sumX1Y_norm - sumX1X2_norm * sumX2Y_norm) / &
   & (sumX1X1_norm * sumX2X2_norm - sumX1X2_norm**2)
b2 = (sumX1X1_norm * sumX2Y_norm - sumX1X2_norm * sumX1Y_norm) / &
   & (sumX1X1_norm * sumX2X2_norm - sumX1X2_norm**2)
a  = (sumY - b1 * sumX1 - b2 * sumX2) / n
write(*,*) b1, b2, a

! Calculate sum of squares, R values
ssReg = b1 * sumX1Y_norm + b2 * sumX2Y_norm
ssRes = sumYY_norm - ssReg
rSquared = ssReg / sumYY_norm
fTest = (ssReg / k) / (ssRes / (n - k - 1))
write(*,*) ssReg, ssRes, rSquared, fTest


! Write out pertinent info to file.
open(96, file = 'regression_output.txt')
   do i = 1, n
      write(96, '(F6.4, F6.4)') y(i), a + b1 * x1(i) + b2 * x2(i)
   end do
close(96)

end program