clear	
clear matrix
set matsize 800
set mem 500m
 
cd ***ADD LOCATION OF DATAFILE***

use ***.csv file as stata .dta file


/// add pitcher & pos player WARs to determine team WARs
	
	gen war=pwar+nonpwar

    // determine win pct of teams
	gen wpct=W/(W+L)


	/// determine WAR-per-game (necessary to reflect varying season-lengths)
		
	gen warpg =war/GP


    /// regress winning pct on War-per-game



	regress wpct warpg


    /// use regression model to determine estimates for winning percentages of teams after computing WAR-per-game by hand
    /// replace the values used in the Margins command with those
	


	margins, at(warpg=(.2864 .2191))
	
// Automatically set the team-strength parameters
local T1_wp = r(table)[1,1]    // First predicted marginal value
local T1_wp_se = r(table)[2,1] // First standard error
local T2_wp = r(table)[1,2]    // Second predicted marginal value
local T2_wp_se = r(table)[2,2] // Second standard error





// Set seed for reproducibility and the number of simulations
set seed 12345
local n_trials = 2500

// Create dataset to store outcomes
clear
set obs `n_trials'
gen T1_series_wins = 0
gen T2_series_wins = 0
gen games_played = 0

// variable to keep track of simulated Team 1 one-game WP
gen T1_1_game_wp=0

// Add variables to track the number of games for each series
gen T1_won_in_4 = 0
gen T1_won_in_5 = 0
gen T1_won_in_6 = 0
gen T1_won_in_7 = 0
gen T2_won_in_4 = 0
gen T2_won_in_5 = 0
gen T2_won_in_6 = 0
gen T2_won_in_7 = 0

// Monte Carlo simulation loop
quietly forvalues i = 1/`n_trials' {
    // Initialize game count and series outcome
    local T1_wins = 0
    local T2_wins = 0
    local games_played = 0

    while (`T1_wins' < 4 & `T2_wins' < 4) {
        // Draw randomly from distibution associated with regression measure of team strengths 
        local T1_wp_game = rnormal(`T1_wp', `T1_wp_se')
        local T2_wp_game = rnormal(`T2_wp', `T2_wp_se')

        // Calculate the one-game probability of Team 1 winning
       local T1_win_game_prob =   (`T1_wp_game'/(1-`T1_wp_game'))/((`T2_wp_game'/(1-`T2_wp_game')))/((`T1_wp_game'/(1-`T1_wp_game'))/((`T2_wp_game'/(1-`T2_wp_game')))+1)

	   //store T1 one-game win probability for first game of series to enable examination of simulated range of values
	    if (`T1_wins' == 0) replace T1_1_game_wp = `T1_win_game_prob' in `i'
     
       
        // Determine the outcome of the game based on the calculated probability
        local result = runiform() < `T1_win_game_prob'
		
		
	 
        if (`result') {
            local T1_wins = `T1_wins' + 1
        }
        else {
            local T2_wins = `T2_wins' + 1
        }

        local games_played = `games_played' + 1
    }

    // Record the result for the trial
    replace T1_series_wins = 1 if `T1_wins' == 4 in `i'
    replace T2_series_wins = 1 if `T2_wins' == 4 in `i'
    replace games_played = `games_played' in `i'
	

    // Track the number of games for each series
    if (`T1_wins' == 4) {
        if (`games_played' == 4) replace T1_won_in_4 = 1 in `i'
        if (`games_played' == 5) replace T1_won_in_5 = 1 in `i'
        if (`games_played' == 6) replace T1_won_in_6 = 1 in `i'
        if (`games_played' == 7) replace T1_won_in_7 = 1 in `i'
    }
    if (`T2_wins' == 4) {
        if (`games_played' == 4) replace T2_won_in_4 = 1 in `i'
        if (`games_played' == 5) replace T2_won_in_5 = 1 in `i'
        if (`games_played' == 6) replace T2_won_in_6 = 1 in `i'
        if (`games_played' == 7) replace T2_won_in_7 = 1 in `i'
    }
}

/// report simulated range of values for Team 1 1-game win probability

sum T1_1_game_wp,d

// report series results
sum T1_series_wins T2_series_wins

// report series lengths
sum T1_won_in_4 T1_won_in_5 T1_won_in_6 T1_won_in_7
sum T2_won_in_4 T2_won_in_5 T2_won_in_6 T2_won_in_7