Portfolio_Simulation_Modeling
Portfolio Simulation Modeling¶
Copyright By PowCoder代写 加微信 powcoder
Download Data¶
Download S&P 500 index data from Yahoo Finance¶
import pandas as pd
import matplotlib.pyplot as plt
import datetime as datetime
from numpy import *
%matplotlib inline
import pandas_datareader as web
!pip install pandas_datareader
import pandas_datareader as web
Set start date, end date and data source (‘Yahoo Finance’, ‘Quandl’, etc.). Download S&P 500 index data from Yahoo Finance.
start_date = datetime.date(1991,1,1)
end_date = datetime.date(2022,1,1)
# Download S&P 500 index data
SnP500_Ddata = web.DataReader(‘^GSPC’,’yahoo’,start_date,end_date)
SnP500_Ddata.head()
High Low Open Close Volume Adj Close
1991-01-02 330.750000 326.450012 330.200012 326.450012 126280000 326.450012
1991-01-03 326.529999 321.899994 326.459991 321.910004 141450000 321.910004
1991-01-04 322.350006 318.869995 321.910004 321.000000 140820000 321.000000
1991-01-07 320.970001 315.440002 320.970001 315.440002 130610000 315.440002
1991-01-08 316.970001 313.790009 315.440002 314.899994 143390000 314.899994
Transform daily data into annual data:¶
# Create a time-series of annual data points from daily data
SnP500_Adata = SnP500_Ddata.resample(‘A’).last()
SnP500_Adata[[‘Volume’,’Adj Close’]].tail()
Volume Adj Close
2017-12-31 2443490000 2673.610107
2018-12-31 3442870000 2506.850098
2019-12-31 2893810000 3230.780029
2020-12-31 3172510000 3756.070068
2021-12-31 2446190000 4766.180176
Compute annual return of S&P 500 index:¶
SnP500_Adata[[‘Adj Close’]] = SnP500_Adata[[‘Adj Close’]].apply(pd.to_numeric, errors=’ignore’)
SnP500_Adata[‘returns’] = SnP500_Adata[‘Adj Close’] / SnP500_Adata[‘Adj Close’].shift(1) -1
SnP500_Adata = SnP500_Adata.dropna()
print (SnP500_Adata[‘returns’])
1992-12-31 0.044643
1993-12-31 0.070552
1994-12-31 -0.015393
1995-12-31 0.341107
1996-12-31 0.202637
1997-12-31 0.310082
1998-12-31 0.266686
1999-12-31 0.195260
2000-12-31 -0.101392
2001-12-31 -0.130427
2002-12-31 -0.233660
2003-12-31 0.263804
2004-12-31 0.089935
2005-12-31 0.030010
2006-12-31 0.136194
2007-12-31 0.035296
2008-12-31 -0.384858
2009-12-31 0.234542
2010-12-31 0.127827
2011-12-31 -0.000032
2012-12-31 0.134057
2013-12-31 0.296012
2014-12-31 0.113906
2015-12-31 -0.007266
2016-12-31 0.095350
2017-12-31 0.194200
2018-12-31 -0.062373
2019-12-31 0.288781
2020-12-31 0.162589
2021-12-31 0.268927
Freq: A-DEC, Name: returns, dtype: float64
Compute average annual return and standard deviation of return for S&P 500 index:¶
SnP500_mean_ret = float(SnP500_Adata[[‘returns’]].mean())
SnP500_std_ret = float(SnP500_Adata[[‘returns’]].std())
print (“S&P 500 average return = %g%%, st. dev = %g%%” % (100*SnP500_mean_ret, 100*SnP500_std_ret))
S&P 500 average return = 9.88999%, st. dev = 16.8951%
Simulation Example 1¶
We want to invest $\$$1000 in the US stock market for 1 year: $v_0 = 1000$
v0 = 1000 # Initial capital
In our example we assume that the return of the market over the next year follow Normal distribution.
Between 1992 and 2022, S&P 500 returned 9.89% per year on average with a standard deviation of 16.90%.
Generate 100 scenarios for the market return over the next year (draw 100 random numbers from a Normal distribution with mean 9.89% and standard deviation of 16.90%):
Ns = 100 # Number of scenarios
r01 = random.normal(SnP500_mean_ret, SnP500_std_ret, Ns)
array([ 0.1375027 , -0.09319938, 0.13307919, 0.35499694, -0.04877715,
0.27564402, 0.11619371, -0.15878737, 0.31654595, 0.03686162,
-0.13288567, 0.28733357, 0.12547818, -0.0276825 , 0.19340088,
0.1677725 , 0.24199379, 0.03088336, 0.10895233, -0.14099148,
0.1977495 , -0.09640502, 0.27998324, 0.28231808, 0.13457243,
-0.2014965 , 0.05497369, 0.17722856, 0.18276941, 0.16527064,
0.40547678, 0.35433638, 0.08711762, 0.03006333, 0.00184905,
0.06033029, -0.0046672 , 0.31979466, 0.0164282 , 0.14794374,
0.02742182, 0.44723226, -0.07484246, 0.3562332 , -0.0291925 ,
-0.14579422, 0.19648547, 0.23792232, 0.25218973, 0.26888589,
-0.05642739, 0.37158784, -0.24566678, 0.09429163, 0.01539461,
0.27825809, 0.04750595, 0.09054314, 0.24449008, 0.33204064,
0.03050717, 0.16483235, -0.07917764, -0.12333402, -0.02137019,
-0.01807691, -0.30444064, 0.13240615, 0.41874921, 0.01823566,
0.01243735, -0.01974399, -0.04185282, -0.07312779, 0.03731275,
0.13205751, -0.04922398, 0.16478846, 0.32611105, 0.14127037,
0.26148578, 0.40534446, 0.29532558, 0.23746478, -0.10464081,
-0.15748783, 0.03942849, 0.20604163, 0.60244983, 0.2590744 ,
-0.24749653, 0.17399643, 0.14862128, -0.11122591, -0.285295 ,
0.01168187, 0.23139099, -0.22145477, -0.13503132, 0.39032693])
Value of investment at the end of year 1:
$v_1 = v_0 + r_{0,1}\cdot v_0 = (1 + r_{0,1})\cdot v_0$
# Distribution of value at the end of year 1
v1 = (r01 + 1) * v0
array([1137.50270418, 906.80062161, 1133.07918539, 1354.99694156,
951.22285108, 1275.64401792, 1116.19370588, 841.21263095,
1316.54595403, 1036.86161635, 867.11432884, 1287.33356827,
1125.47817856, 972.31749904, 1193.40087712, 1167.77249639,
1241.99378815, 1030.88335965, 1108.95232776, 859.00852039,
1197.74949624, 903.5949815 , 1279.98324197, 1282.31807964,
1134.57243231, 798.50349914, 1054.97368732, 1177.22855615,
1182.76940926, 1165.27064342, 1405.47677579, 1354.33638 ,
1087.11762217, 1030.06333334, 1001.84905404, 1060.33029443,
995.33280177, 1319.79466346, 1016.42819501, 1147.94373965,
1027.42181585, 1447.23226171, 925.15754183, 1356.23319938,
970.80750257, 854.20577635, 1196.485465 , 1237.92232191,
1252.18972769, 1268.88589104, 943.57261202, 1371.58783516,
754.33321671, 1094.2916341 , 1015.39461225, 1278.25808982,
1047.50594755, 1090.54314192, 1244.49008256, 1332.04064207,
1030.50717454, 1164.83234733, 920.82235806, 876.66597732,
978.62981076, 981.92308596, 695.55935719, 1132.40615281,
1418.74920544, 1018.23565512, 1012.43735065, 980.25601034,
958.14718304, 926.87221184, 1037.31274612, 1132.05750707,
950.77602158, 1164.78845544, 1326.11104846, 1141.27036606,
1261.48577812, 1405.34446228, 1295.32557691, 1237.46477843,
895.35918966, 842.51216511, 1039.42848799, 1206.04162938,
1602.44983468, 1259.07439957, 752.5034667 , 1173.99643296,
1148.62128371, 888.77409397, 714.70500223, 1011.68186759,
1231.39098568, 778.54522605, 864.96868338, 1390.32693088])
Compute statistical measures from the distribution¶
1100.7687167764386
Standard deviation:
183.72653726457224
Minimum, maximum:
min(v1), max(v1)
(695.5593571872389, 1602.4498346760681)
Compute percentiles/quantiles¶
5th percentile, median, 95th percentile:
percentile(v1, [5, 50,95])
array([ 797.50558549, 1112.57301682, 1391.07780745])
sortedScen = sorted(v1) # sort scenarios
mean(sortedScen[4:6]), mean(sortedScen[49:51]), mean(sortedScen[94:96])
(788.5243625973617, 1112.573016819433, 1397.8356965781786)
Alternative way to compute percentiles/quantiles
sortedScen[int(Ns-(1-0.05)*Ns)-1] # 5th percentile
778.5452260516112
# Alternative way to compute percentiles/quantiles
sortedScen[int(Ns-(1-0.95)*Ns)-1] # 95th percentile
1390.3269308750187
Visualize results¶
Plot a histogram of the distribution of outcomes for v1:
hist, bins = histogram(v1)
positions = (bins[:-1] + bins[1:]) / 2
plt.bar(positions, hist, width=60)
plt.xlabel(‘portfolio value in 1 year’)
plt.ylabel(‘frequency’)
plt.show()
Simulated paths over time:
# Plot simulated paths over time
for res in v1:
plt.plot((0,1), (v0, res))
plt.xlabel(‘time step’)
plt.ylabel(‘portfolio value’)
plt.title(‘Simulated Value Paths’)
plt.show()
Simulation Example 2¶
You are planning for retirement and decide to invest in the market for the next 30 years (instead of only the next year as in example 1).
Assume that every year your investment returns from investing into the S&P 500 will follow a Normal distribution with the mean and standard deviation as in example 1.
Your initial capital is still \$1000
v0 = 1000 # Initial capital
Between 1992 and 2022, S&P 500 returned 9.89% per year on average with a standard deviation of 16.90%.
Simulate 30 columns of 100 observations each of single period returns:
r_speriod30 = random.normal(SnP500_mean_ret, SnP500_std_ret, (Ns, 30))
r_speriod30
array([[ 0.12833353, -0.06821577, -0.30048505, …, -0.08739151,
0.11830901, 0.16103138],
[ 0.14956169, 0.37917207, 0.10713295, …, 0.06578568,
0.16694179, 0.42415946],
[-0.00161962, 0.15484477, 0.25853146, …, 0.11377973,
0.4734983 , 0.0077709 ],
[-0.03108745, 0.09881029, 0.169197 , …, 0.03291492,
0.23706931, -0.04028601],
[ 0.11647758, 0.13916284, -0.00413712, …, -0.06907805,
0.03893452, 0.15268992],
[ 0.04939218, 0.28152746, -0.01211741, …, 0.11857633,
0.06563149, 0.06619776]])
Compute and plot $v_{30}$
v30 = prod(1 + r_speriod30 , 1) * v0
hist, bins = histogram(v30)
positions = (bins[:-1] + bins[1:]) / 2
width = (bins[1] – bins[0]) * 0.8
plt.bar(positions, hist, width=width)
plt.xlabel(‘portfolio value after 30 years’)
plt.ylabel(‘frequency’)
plt.show()
Simulated paths over time:
for scenario in r_speriod30:
y = [prod(1 + scenario[0:i]) * v0 for i in range(0,31)]
plt.plot(range(0,31), y)
plt.xlabel(‘time step’)
plt.ylabel(‘portfolio value’)
plt.title(‘Simulated Value Paths’)
plt.show()
Compute percentiles/quantiles¶
5th percentile, median, 95th percentile:
percentile(v30, [5, 50,95])
array([ 2578.85977252, 10914.23236257, 40075.24046231])
sortedScen = sorted(v30) # sort scenarios
sortedScen[int(Ns-(1-0.05)*Ns)-1], sortedScen[int(Ns-(1-0.95)*Ns)-1]
(2138.2941333117014, 39905.185316077426)
Simulation Example 3¶
Download US Treasury bill data from Federal Reserve:¶
# Download 3-month T-bill rates from Federal Reserve
start_date_b = datetime.date(1992,1,1)
end_date_b = datetime.date(2022,1,1)
TBill_Ddata = web.DataReader(‘DTB3′,’fred’,start_date_b,end_date_b)
TBill_Ddata.head()
1992-01-01 NaN
1992-01-02 3.86
1992-01-03 3.85
1992-01-06 3.81
1992-01-07 3.75
Transform daily data into annual data:¶
# Create a time-series of annual data points from daily data
TBill_Adata = TBill_Ddata.resample(‘A’).last()
TBill_Adata[[‘DTB3’]].tail()
2017-12-31 1.37
2018-12-31 2.40
2019-12-31 1.52
2020-12-31 0.09
2021-12-31 0.06
Compute annual return for bonds:¶
TBill_Adata[‘returns’] = TBill_Adata[‘DTB3’] / 100
TBill_Adata = TBill_Adata.dropna()
print (TBill_Adata[‘returns’])
1992-12-31 0.0308
1993-12-31 0.0301
1994-12-31 0.0553
1995-12-31 0.0496
1996-12-31 0.0507
1997-12-31 0.0522
1998-12-31 0.0437
1999-12-31 0.0517
2000-12-31 0.0573
2001-12-31 0.0171
2002-12-31 0.0120
2003-12-31 0.0093
2004-12-31 0.0218
2005-12-31 0.0399
2006-12-31 0.0489
2007-12-31 0.0329
2008-12-31 0.0011
2009-12-31 0.0006
2010-12-31 0.0012
2011-12-31 0.0002
2012-12-31 0.0005
2013-12-31 0.0007
2014-12-31 0.0004
2015-12-31 0.0016
2016-12-31 0.0050
2017-12-31 0.0137
2018-12-31 0.0240
2019-12-31 0.0152
2020-12-31 0.0009
2021-12-31 0.0006
Freq: A-DEC, Name: returns, dtype: float64
Compute average annual return and standard deviation of return for bonds:¶
TBill_mean_ret = float(TBill_Adata[[‘returns’]].mean())
TBill_std_ret = float(TBill_Adata[[‘returns’]].std())
print (“T-bill average return = %g%%, st. dev = %g%%” % (100*TBill_mean_ret, 100*TBill_std_ret))
T-bill average return = 2.23%, st. dev = 2.09322%
Compute covariance matrix:¶
covMat = cov(array(SnP500_Adata[[‘returns’]]),array(TBill_Adata[[‘returns’]]),rowvar=0)
array([[0.0285443 , 0.00036953],
[0.00036953, 0.00043816]])
Simulate portfolio:¶
v0 = 1000 # Initial capital
Ns = 5000 # Number of scenarios
mu = [SnP500_mean_ret, TBill_mean_ret] # Expected return
[0.09889988936460865, 0.022300000000000004]
stockRet = ones(Ns)
bondsRet = ones(Ns)
scenarios = random.multivariate_normal(mu, covMat, Ns)
for year in range(1, 31):
scenarios = random.multivariate_normal(mu, covMat, Ns)
stockRet *= (1 + scenarios[:,0])
bondsRet *= (1 + scenarios[:,1])
v30 = 0.5 * v0 * stockRet + 0.5 * v0 * bondsRet
hist, bins = histogram(v30, bins = 100)
positions = (bins[:-1] + bins[1:]) / 2
width = (bins[1] – bins[0]) * 0.8
plt.bar(positions, hist, width=width)
plt.xlabel(‘portfolio value after 30 years’)
plt.ylabel(‘frequency’)
plt.show()
Simulation Example 4¶
Compare two portfolios
# Compute portfolios by iterating through different combinations of weights
v30comp = []
for w in arange(0.2, 1.01, 0.2):
v30comp += [w * v0 * stockRet + (1 – w) * v0 * bondsRet]
# Plot a histogram of the distribution of
# differences in outcomes for v30
# (Stratery 4 – Strategy 2)
v30d = v30comp[3] – v30comp[1]
hist, bins = histogram(v30d, bins = 50)
positions = (bins[:-1] + bins[1:]) / 2
width = (bins[1] – bins[0]) * 0.8
plt.bar(positions, hist, width=width)
plt.show()
# Compute number of elements in v30d that are > 0 and < 0 and compare
pos_count = (v30d > 0).sum()
neg_count = (v30d <= 0).sum()
print (u"""Strategy 4 was better in %d cases (%d%%)
Strategy 2 was better in %d cases (%d%%)
Difference = %d scenarios""" % (pos_count, round(100*(pos_count/Ns)), neg_count, round(100*(neg_count/Ns)), pos_count - neg_count))
Strategy 4 was better in 4882 cases (98%)
Strategy 2 was better in 118 cases (2%)
Difference = 4764 scenarios
程序代写 CS代考 加微信: powcoder QQ: 1823890830 Email: powcoder@163.com