!  OLS.f90 
INCLUDE 'link_f90_static.h'

module globals

    implicit none
    save
    
    integer, parameter:: dble = kind(1d0)
    
    integer, parameter:: N = 100
    real(dble), parameter:: alpha = 5.0d0
    real(dble), parameter:: beta = 1.0d0
    
    real(dble) xvect(N), yvect(N)
    
end module globals

module subroutines
    use globals
    use lin_sol_gen_int
    use rnnoa_int
    use neqnf_int
    use umpol_int

    contains
    
    subroutine GenerateData(alpha, beta, N)
        implicit none
        
        real(dble), intent(in):: alpha, beta
        integer, intent(in):: N
        
        real(dble) epsVect(N)
        integer i
        
        xvect = (/(i, i=1,N)/)
        call rnnoa(epsVect)  !Generate standard normal random numbers
        yvect = alpha + beta * xvect + epsVect
        
    end subroutine GenerateData
    
    subroutine OLS1(alphaHat, betaHat)
        implicit none
        
        real(dble), intent(out):: alphaHat, betaHat
        
        real(dble) Xmat(N,2), bvect(2,1), XpX(2,2), Xpy(2,1)
        integer i
        
        Xmat(:,1) = 1.0d0
        Xmat(:,2) = xvect
        XpX = matmul(transpose(Xmat),Xmat)
        Xpy(:,1) = matmul(transpose(Xmat), yvect)
      
        call lin_sol_gen(XpX, Xpy, bvect)
        
        alphaHat = bvect(1,1)
        betaHat = bvect(2,1)
    
    end subroutine OLS1

    subroutine OLS2(alphaHat, betaHat)
        implicit none
        
        real(dble), intent(out):: alphaHat, betaHat
        
        real(dble) meanX, meanY

        meanX = sum(xvect)/N
        meanY = sum(yvect)/N

        betaHat = sum((xvect-meanX)*(yvect-meanY))/sum((xvect-meanX)**2)
        alphaHat = meanY - betaHat * meanX      

    end subroutine OLS2
    
    subroutine OLS3(alphaHat, betaHat)
        implicit none
        
        real(dble), intent(out):: alphaHat, betaHat
        
        real(dble) sol(2), guess(2), fnorm
        real(dble):: errrel = 1.e-10
        external focs
        
        guess = (/ 1.0d0, 2.0d0 /)
        
        call neqnf(focs, sol, xguess = guess, fnorm = fnorm, errrel = errrel) 
        !print '(a,es10.2)', "fnorm =", fnorm
        
        alphaHat = sol(1)
        betaHat = sol(2)
                
    end subroutine OLS3
    
    subroutine OLS4(alphaHat, betaHat)
        implicit none
        
        real(dble), intent(out):: alphaHat, betaHat
        
        real(dble) sol(2), guess(2)
        real(dble):: ftol = 1.e-10
        external sse
        
        guess = (/ 1.0d0, 2.0d0 /)
        
        call umpol(sse, sol, xguess = guess, ftol = ftol) 
        
        alphaHat = sol(1)
        betaHat = sol(2)
                
    end subroutine OLS4    
    
end module subroutines

subroutine focs(in, out, M)
    use globals
    implicit none
    real(dble) in(M), out(M)
    integer M
    
    real(dble) alphaHat, betaHat
    
    alphaHat = in(1)
    betaHat = in(2)
    
    out(1) = sum(yvect-alphaHat-betaHat*xvect)
    out(2) = sum((yvect-alphaHat-betaHat*xvect)*xvect)
    
end subroutine focs  

subroutine sse(M, in, out)
    use globals
    implicit none
    real(dble) in(M), out
    integer M
    
    real(dble) alphaHat, betaHat
    
    alphaHat = in(1)
    betaHat = in(2)
    
    out = sum( ( yvect- alphaHat - betaHat * xvect )**2 )
    
end subroutine sse  

program OLSExample

    use subroutines

    implicit none
    
    real(dble) alphaHat, betaHat
    integer i, j

    call GenerateData(alpha, beta, N)
    
    call OLS1(alphaHat, betaHat)
    print 101, "OLS1: alphaHat =", alphaHat, " betaHat =", betaHat 
    
    call OLS2(alphaHat, betaHat)
    print 101, "OLS2: alphaHat =", alphaHat, " betaHat =", betaHat 
    
    call OLS3(alphaHat, betaHat)
    print 101, "OLS3: alphaHat =", alphaHat, " betaHat =", betaHat 
    
    call OLS4(alphaHat, betaHat)
    print 101, "OLS4: alphaHat =", alphaHat, " betaHat =", betaHat     

    101 format (a,f7.4,a,f7.4) 
end program OLSExample