Program financial_risk_analysis
    Use, Intrinsic :: iso_fortran_env, Only: wp => real64
    Use :: nag_library, Only: g05kff, g05skf
    Use :: timer_module
    Implicit None

    ! パラメータ
    Integer, Parameter :: num_assets = 1000
    Integer, Parameter :: time_horizon = 252
    Integer, Parameter :: total_simulations = 1000000
    Integer, Parameter :: target_batch_size = 2000

    ! 変数
    Real (wp), Allocatable :: returns(:, :), weights(:), portfolio_values(:)
    Real (wp), Allocatable :: random_numbers(:)
    Real (wp) :: mu, sigma, initial_value, var, es
    Integer :: ifail
    Integer, Allocatable :: state(:)
    Integer :: num_batches, batch_size, simulations_per_batch

    Call init_timer()

    ! バッチサイズの計算
    simulations_per_batch = min(target_batch_size, total_simulations)
    batch_size = simulations_per_batch*num_assets*time_horizon

    ! 乱数生成器の初期化
    Call initialize_random_generator(.False.) ! .True.なら再現可能乱数、.False.なら毎回違う乱数

    ! データの初期化
    Call initialize_data()
    Call elapsed_from_last_check('Initialization Completed')

    ! シミュレーションの実行
    Call run_simulations()

    Call elapsed_from_last_check('Simulation Completed')

    ! リスク指標の計算
    Call calculate_risk_metrics()

    Call elapsed_from_last_check('Metrics Calculation Completed')

    ! 結果の出力
    Call output_results()

    Call elapsed_total('Total Execution Time')

Contains

    Subroutine initialize_random_generator(repeatable)
        Logical, Intent (In) :: repeatable
        Integer :: genid, subid, lseed, lstate, ifail
        Integer, Allocatable :: seed(:)

        genid = 3                  ! Mersenne Twister
        subid = 0
        lstate = 0
        ifail = 0

        Allocate (state(lstate))

        If (repeatable) Then
            lseed = 1
            Allocate (seed(lseed))
            seed(1) = 123456       ! 固定シードを設定

            ! 最初の呼び出しで必要なlstateのサイズを取得
            Call g05kff(genid, subid, seed, lseed, state, lstate, ifail)

            ! stateを割り当て
            If (allocated(state)) Deallocate (state)
            Allocate (state(lstate))

            ! 実際の初期化
            Call g05kff(genid, subid, seed, lseed, state, lstate, ifail)

            Print *, 'Repeatable random generator initialized with seed:', seed(1)
            Deallocate (seed)
        Else
            ! 非リピータブルな初期化
            Call g05kgf(genid, subid, state, lstate, ifail)

            ! stateを割り当て
            If (allocated(state)) Deallocate (state)
            Allocate (state(lstate))

            ! 実際の初期化
            Call g05kgf(genid, subid, state, lstate, ifail)

            Print *, 'Non-repeatable random generator initialized.'
        End If
    End Subroutine initialize_random_generator


    Subroutine initialize_data()
        Integer :: ifail, n_total
        Real (wp) :: mu_returns, sigma_returns
        Integer :: genid, subid, lseed, lstate
        Integer, Allocatable :: seed(:), state_returns(:)

        Allocate (returns(num_assets,time_horizon))
        Allocate (weights(num_assets))
        Allocate (portfolio_values(total_simulations))
        Allocate (random_numbers(batch_size))

        ! Parameters for the normal distribution of returns
        mu_returns = 0.0_wp        ! Mean return
        sigma_returns = 0.01_wp    ! Standard deviation set to 1%

        n_total = num_assets*time_horizon
        ifail = 0

        genid = 3                  ! Mersenne Twister
        subid = 0
        lstate = 0
        ifail = 0

        Allocate (state_returns(lstate))
        lseed = 1
        Allocate (seed(lseed))
        seed(1) = 54321           ! returnsは固定シードで生成

        ! 最初の呼び出しで必要なlstateのサイズを取得
        Call g05kff(genid, subid, seed, lseed, state_returns, lstate, ifail)
        If (allocated(state_returns)) Deallocate (state_returns)
        Allocate (state_returns(lstate))
        ! 実際の初期化
        Call g05kff(genid, subid, seed, lseed, state_returns, lstate, ifail)
        Deallocate (seed)


        ! Generate normally distributed returns for assets using state_returns
        Call g05skf(n_total, mu_returns, sigma_returns**2, state_returns, returns, ifail)

        ! Deallocate state_returns
        Deallocate (state_returns)

        ! Initialize weights: equal-weighted portfolio
        weights = 1.0_wp/real(num_assets, wp)

        ! Initial portfolio value
        initial_value = 1000000.0_wp

        ! Parameters for the random shocks added in simulations
        mu = 0.0_wp
        sigma = 0.01_wp            ! Standard deviation remains 1%

        Print *, 'Data initialization completed with normal distribution for returns.'
        Print *, 'First few returns:', returns(1:3, 1)
        Print *, 'Sum of weights:', sum(weights)
        Print *, 'Initial portfolio value:', initial_value
        Print *, 'Mean of returns:', sum(returns)/size(returns)
        Print *, 'Std dev of returns:', sqrt(sum((returns-sum(returns)/size(returns))**2)/(size(returns)-1))
    End Subroutine initialize_data


    Subroutine run_simulations()
        Integer :: batch

        ! バッチサイズのチェック
        If (mod(total_simulations,simulations_per_batch)/=0) Then
            Print *, 'Error: Total simulations must be divisible by simulations per batch.'
            Print *, 'Total simulations:', total_simulations
            Print *, 'Simulations per batch:', simulations_per_batch
            Stop 1                 ! プログラムを終了し、エラーコード1を返す
        End If

        num_batches = total_simulations/simulations_per_batch

        Do batch = 1, num_batches
            Call g05skf(int(batch_size), mu, sigma**2, state, random_numbers, ifail)
            Call process_batch(random_numbers, batch)
            Print '(A,I0,A,I0,A,F6.2,A)', 'Batch ', batch, ' of ', num_batches, ' completed (', 100.0*real(batch)/real(num_batches), &
                '%)'
        End Do

    End Subroutine run_simulations


    Subroutine process_batch(random_nums, batch_num)
        Real (wp), Intent (In) :: random_nums(:)
        Integer, Intent (In) :: batch_num
        Integer :: sim, asset, time, random_index
        Real (wp) :: portfolio_value, portfolio_return, asset_return
        Integer :: start_sim, end_sim

        start_sim = (batch_num-1)*simulations_per_batch + 1
        end_sim = min(start_sim+simulations_per_batch-1, total_simulations)

        Do sim = start_sim, end_sim
            portfolio_value = initial_value
            random_index = (sim-start_sim)*num_assets*time_horizon + 1

            Do time = 1, time_horizon
                portfolio_return = 0.0_wp
                Do asset = 1, num_assets
                    ! Corrected: Remove multiplication by sigma
                    asset_return = returns(asset, time) + random_nums(random_index)
                    random_index = random_index + 1

                    ! Add to portfolio return
                    portfolio_return = portfolio_return + weights(asset)*asset_return
                End Do

                ! Update portfolio value
                portfolio_value = portfolio_value*(1.0_wp+portfolio_return)

            End Do

            portfolio_values(sim) = portfolio_value

        End Do
    End Subroutine process_batch


    Subroutine calculate_risk_metrics()
        Integer :: i
        Real (wp) :: confidence_level = 0.95_wp
        Integer :: var_index
        Real (wp) :: loss, total_loss, average_loss

        ! Sort the portfolio values in ascending order
        Call sort_ascending(portfolio_values)

        ! Calculate the index for VaR
        var_index = floor((1.0_wp-confidence_level)*total_simulations) + 1

        ! Ensure the index does not exceed array bounds
        If (var_index>total_simulations) var_index = total_simulations

        ! Calculate VaR
        var = max(initial_value-portfolio_values(var_index), 0.0_wp)

        ! Calculate ES
        total_loss = 0.0_wp
        Do i = 1, var_index
            loss = max(initial_value-portfolio_values(i), 0.0_wp)
            total_loss = total_loss + loss
        End Do
        average_loss = total_loss/real(var_index, wp)
        es = average_loss

    End Subroutine calculate_risk_metrics


    Recursive Subroutine sort_ascending(arr)
        Real (wp), Intent (Inout) :: arr(:)
        Integer :: i, j
        Real (wp) :: pivot, temp

        If (size(arr)>1) Then
            pivot = arr(size(arr)/2)
            i = 1
            j = size(arr)
            Do
                Do While (arr(i)<pivot)
                    i = i + 1
                End Do
                Do While (arr(j)>pivot)
                    j = j - 1
                End Do
                If (i>=j) Exit
                temp = arr(i)
                arr(i) = arr(j)
                arr(j) = temp
                i = i + 1
                j = j - 1
            End Do
            Call sort_ascending(arr(:j))
            Call sort_ascending(arr(j+1:))
        End If
    End Subroutine sort_ascending


    Subroutine output_results()
        Real (wp) :: var_loss, var_percentage, es_loss, es_percentage

        ! Corrected: Use var and es directly
        var_loss = var
        var_percentage = (var_loss/initial_value)*100.0_wp
        es_loss = es
        es_percentage = (es_loss/initial_value)*100.0_wp

        Print '(A,F12.2,A)', 'Initial Portfolio Value: ', initial_value, ' yen'
        Print '(A,F12.2,A,F5.2,A)', 'Value at Risk (VaR): ', var_loss, ' yen (', var_percentage, '% of initial investment)'
        Print '(A,F12.2,A,F5.2,A)', 'Expected Shortfall (ES): ', es_loss, ' yen (', es_percentage, '% of initial investment)'
        Print *, 'Time Horizon: 252 days (1 year)'
        Print *, 'Confidence Level: 95%'
    End Subroutine output_results


End Program financial_risk_analysis