Simulation and Modeling: Building Computational Models with Python

Simulation and modeling enable understanding complex systems through computational experiments. Python provides excellent tools for building and analyzing models.

Monte Carlo Simulations

Basic Monte Carlo Method

import numpy as np
import matplotlib.pyplot as plt

# Estimate π using Monte Carlo
def estimate_pi(num_samples):
    """Estimate π by randomly sampling points in a square"""
    # Generate random points in [0,1] x [0,1]
    x = np.random.uniform(0, 1, num_samples)
    y = np.random.uniform(0, 1, num_samples)
    
    # Calculate distance from origin
    distances = np.sqrt(x**2 + y**2)
    
    # Count points inside unit circle
    inside_circle = np.sum(distances <= 1)
    
    # Estimate π
    pi_estimate = 4 * inside_circle / num_samples
    
    return pi_estimate, x, y, distances

# Run simulation
num_samples = 100000
pi_est, x, y, distances = estimate_pi(num_samples)

print(f"Estimated π: {pi_est:.6f}")
print(f"Actual π: {np.pi:.6f}")
print(f"Error: {abs(pi_est - np.pi):.6f}")

# Visualize
plt.figure(figsize=(10, 8))
inside = distances <= 1
plt.scatter(x[inside], y[inside], c='red', s=1, alpha=0.5, label='Inside circle')
plt.scatter(x[~inside], y[~inside], c='blue', s=1, alpha=0.5, label='Outside circle')
plt.xlim(0, 1)
plt.ylim(0, 1)
plt.gca().set_aspect('equal')
plt.legend()
plt.title(f'Monte Carlo Estimation of π ≈ {pi_est:.4f}')
plt.show()

Portfolio Risk Analysis

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

def simulate_portfolio_returns(initial_investment, annual_return, annual_volatility, years, num_simulations):
    """Simulate portfolio returns using geometric Brownian motion"""
    
    # Time parameters
    dt = 1/252  # Daily time step
    num_steps = years * 252
    
    # Initialize results
    simulations = np.zeros((num_simulations, num_steps))
    simulations[:, 0] = initial_investment
    
    # Run simulations
    for i in range(1, num_steps):
        # Random returns
        random_returns = np.random.normal(
            annual_return * dt,
            annual_volatility * np.sqrt(dt),
            num_simulations
        )
        
        # Update portfolio values
        simulations[:, i] = simulations[:, i-1] * (1 + random_returns)
    
    return simulations

# Parameters
initial_investment = 100000
annual_return = 0.08
annual_volatility = 0.15
years = 10
num_simulations = 1000

# Run simulation
results = simulate_portfolio_returns(
    initial_investment, annual_return, annual_volatility, years, num_simulations
)

# Analyze results
final_values = results[:, -1]
percentiles = np.percentile(final_values, [5, 25, 50, 75, 95])

print(f"Initial investment: ${initial_investment:,.0f}")
print(f"Expected final value: ${np.mean(final_values):,.0f}")
print(f"5th percentile: ${percentiles[0]:,.0f}")
print(f"25th percentile: ${percentiles[1]:,.0f}")
print(f"Median: ${percentiles[2]:,.0f}")
print(f"75th percentile: ${percentiles[3]:,.0f}")
print(f"95th percentile: ${percentiles[4]:,.0f}")

# Visualize
plt.figure(figsize=(12, 6))
plt.plot(results.T, alpha=0.1, color='blue')
plt.plot(np.mean(results, axis=0), color='red', linewidth=2, label='Mean')
plt.xlabel('Days')
plt.ylabel('Portfolio Value ($)')
plt.title('Portfolio Value Simulations')
plt.legend()
plt.grid(True)
plt.show()

Agent-Based Modeling

Simple Agent-Based Model

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

class Agent:
    def __init__(self, x, y, velocity):
        self.x = x
        self.y = y
        self.vx = velocity[0]
        self.vy = velocity[1]
    
    def move(self, width, height):
        """Move agent and handle boundaries"""
        self.x += self.vx
        self.y += self.vy
        
        # Bounce off walls
        if self.x <= 0 or self.x >= width:
            self.vx *= -1
        if self.y <= 0 or self.y >= height:
            self.vy *= -1
        
        # Keep within bounds
        self.x = np.clip(self.x, 0, width)
        self.y = np.clip(self.y, 0, height)
    
    def distance_to(self, other):
        """Calculate distance to another agent"""
        return np.sqrt((self.x - other.x)**2 + (self.y - other.y)**2)

class ABMSimulation:
    def __init__(self, num_agents, width, height):
        self.agents = []
        self.width = width
        self.height = height
        
        # Create agents
        for _ in range(num_agents):
            x = np.random.uniform(0, width)
            y = np.random.uniform(0, height)
            vx = np.random.uniform(-1, 1)
            vy = np.random.uniform(-1, 1)
            self.agents.append(Agent(x, y, (vx, vy)))
    
    def step(self):
        """Execute one simulation step"""
        for agent in self.agents:
            agent.move(self.width, self.height)
    
    def get_positions(self):
        """Get all agent positions"""
        positions = np.array([[a.x, a.y] for a in self.agents])
        return positions

# Run simulation
sim = ABMSimulation(num_agents=50, width=100, height=100)

# Collect data
positions_history = []
for _ in range(100):
    sim.step()
    positions_history.append(sim.get_positions().copy())

# Visualize
fig, ax = plt.subplots(figsize=(8, 8))

def animate(frame):
    ax.clear()
    positions = positions_history[frame]
    ax.scatter(positions[:, 0], positions[:, 1], s=50, alpha=0.6)
    ax.set_xlim(0, 100)
    ax.set_ylim(0, 100)
    ax.set_title(f'Agent-Based Model - Step {frame}')
    ax.set_aspect('equal')

# Create animation
anim = FuncAnimation(fig, animate, frames=len(positions_history), interval=50)
plt.show()

Epidemic Simulation (SIR Model)

import numpy as np
import matplotlib.pyplot as plt

class SIRModel:
    def __init__(self, population, initial_infected, transmission_rate, recovery_rate):
        self.population = population
        self.susceptible = population - initial_infected
        self.infected = initial_infected
        self.recovered = 0
        self.transmission_rate = transmission_rate
        self.recovery_rate = recovery_rate
        
        # History
        self.history = {
            'susceptible': [self.susceptible],
            'infected': [self.infected],
            'recovered': [self.recovered]
        }
    
    def step(self):
        """Execute one time step"""
        # New infections
        contact_rate = self.transmission_rate * self.infected * self.susceptible / self.population
        new_infections = np.random.poisson(contact_rate)
        new_infections = min(new_infections, self.susceptible)
        
        # Recoveries
        new_recoveries = np.random.binomial(self.infected, self.recovery_rate)
        
        # Update compartments
        self.susceptible -= new_infections
        self.infected += new_infections - new_recoveries
        self.recovered += new_recoveries
        
        # Record history
        self.history['susceptible'].append(self.susceptible)
        self.history['infected'].append(self.infected)
        self.history['recovered'].append(self.recovered)
    
    def run(self, num_steps):
        """Run simulation for specified steps"""
        for _ in range(num_steps):
            self.step()
    
    def plot(self):
        """Plot results"""
        plt.figure(figsize=(12, 6))
        plt.plot(self.history['susceptible'], label='Susceptible', linewidth=2)
        plt.plot(self.history['infected'], label='Infected', linewidth=2)
        plt.plot(self.history['recovered'], label='Recovered', linewidth=2)
        plt.xlabel('Time (days)')
        plt.ylabel('Number of people')
        plt.title('SIR Epidemic Model')
        plt.legend()
        plt.grid(True)
        plt.show()

# Run simulation
population = 10000
initial_infected = 10
transmission_rate = 0.5
recovery_rate = 0.1

model = SIRModel(population, initial_infected, transmission_rate, recovery_rate)
model.run(200)
model.plot()

# Print statistics
print(f"Final susceptible: {model.susceptible}")
print(f"Final infected: {model.infected}")
print(f"Final recovered: {model.recovered}")
print(f"Attack rate: {model.recovered / population * 100:.1f}%")

Numerical Simulations

Heat Diffusion

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

def heat_diffusion_2d(width, height, num_steps, dt=0.01, diffusion=0.1):
    """Simulate 2D heat diffusion"""
    
    # Initialize temperature grid
    T = np.zeros((height, width))
    T[height//2, width//2] = 100  # Heat source
    
    # Store history
    history = [T.copy()]
    
    # Simulation
    for step in range(num_steps):
        T_new = T.copy()
        
        # Laplacian (finite difference)
        for i in range(1, height-1):
            for j in range(1, width-1):
                laplacian = (T[i+1, j] + T[i-1, j] + T[i, j+1] + T[i, j-1] - 4*T[i, j])
                T_new[i, j] = T[i, j] + diffusion * dt * laplacian
        
        T = T_new
        history.append(T.copy())
    
    return history

# Run simulation
width, height = 50, 50
history = heat_diffusion_2d(width, height, num_steps=100)

# Visualize
fig, axes = plt.subplots(2, 2, figsize=(10, 10))
steps = [0, 25, 50, 99]

for idx, (ax, step) in enumerate(zip(axes.flat, steps)):
    im = ax.imshow(history[step], cmap='hot')
    ax.set_title(f'Step {step}')
    plt.colorbar(im, ax=ax)

plt.tight_layout()
plt.show()

Best Practices

Validate models: Compare with analytical solutions when possible
Sensitivity analysis: Test how results change with parameter variations
Reproducibility: Set random seeds for reproducible results
Efficiency: Use NumPy vectorization for large simulations
Visualization: Visualize results to understand model behavior
Documentation: Document model assumptions and parameters

Common Pitfalls

Bad Practice:

# Don't: No random seed
results = simulate_model()

# Don't: Inefficient loops
for i in range(n):
    for j in range(m):
        result[i, j] = expensive_calculation()

# Don't: Ignore numerical stability
T_new = T + large_dt * laplacian

Good Practice:

# Do: Set random seed
np.random.seed(42)
results = simulate_model()

# Do: Use vectorization
result = expensive_calculation_vectorized()

# Do: Use appropriate time steps
T_new = T + stable_dt * laplacian

Conclusion

Simulation and modeling enable understanding complex systems through computational experiments. Master Monte Carlo methods, agent-based modeling, and numerical simulations to build powerful computational models. Combine with visualization and analysis to extract insights from simulations.