Chapter6 金融與統(tǒng)計(jì)模型

數(shù)據(jù)分析與可視化第六章金融與統(tǒng)計(jì)模型6.1確定性模型6.2隨機(jī)模型6.3閾值模型第六章金融和統(tǒng)計(jì)模型2《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)6.1.1總回報(bào)模型6.1確定性模型3《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)金融時(shí)間序列：金融市場中常見的價(jià)格數(shù)據(jù)等常見金融時(shí)間序列模型：（1）確定性模型：回報(bào)率模型（2）隨機(jī)模型：模特卡洛模擬（3）波動率模型6.1確定性模型4《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)投資收益取決于價(jià)格變化+持有資產(chǎn)數(shù)量總回報(bào)模型：考慮資產(chǎn)價(jià)格隨時(shí)間變化的模型6.1確定性模型：總回報(bào)5《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)

考慮通脹：（1）何時(shí)資產(chǎn)翻倍？（2）回報(bào)率曲線？（3）每年通脹率如何確定？利用代碼（下頁）： Q：初始投資為1萬美元，回報(bào)率為6%，多少年后投資會翻倍

？6《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)6.1確定性模型：總回報(bào)

例：使用matplotlib和NumPy來進(jìn)行探索7《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)importmatplotlib.pyplotaspltprinciple_value=10000#investedamountgrossReturn=1.06#Rtreturn_amt=[]x=[]y=[10000]year=2010return_amt.append(principle_value)x.append(year)6.1確定性模型：總回報(bào)8《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)foriinrange(1,15):return_amt.append(return_amt[i-1]*grossReturn)print("Year-",i,"Returned:",return_amt[i])year+=1x.append(year)y.append(833.33*(year-2010)+principle_value)#setthegridtoappearplt.grid()#plotthereturnvaluescurveplt.plot(x,return_amt,color='r')plt.plot(x,y,color='b')plt.show()6.1確定性模型：總回報(bào)9《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)繪制投資回報(bào)曲線投資收益圖

6.1確定性模型：總回報(bào)例：抵押貸款，35萬美元的貸款金額，利率為5%，期限為30年

10《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)6.1確定性模型：總回報(bào)fromdecimalimportDecimalimportmatplotlib.pyplotasplt

colors=[(31,119,180),(174,199,232),(255,128,0),(255,15,14),(44,160,44),(152,223,138),(214,39,40),(255,173,61),(148,103,189),(197,176,213),(140,86,75),(196,156,148),(227,119,194),(247,182,210),(127,127,127),(199,199,199),(188,189,34),(219,219,141),(23,190,207),(158,218,229)]#ScaletheRGBvaluestothe[0,1]range,whichistheformatmatplotlibaccepts.foriinrange(len(colors)):r,g,b=colors[i]colors[i]=(r/255.,g/255.,b/255.)11《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)6.1確定性模型：總回報(bào)defprintHeaders(term,extra):#printheadersprint("\nExtra-Payment:$"+str(extra)+"Term:"+str(term)+"years")print("---------------------------------------------------------")print('Pmtno'.rjust(6),'','Beg.bal.'.ljust(13),'',)print('Payment'.ljust(9),'','Principal'.ljust(9),'',)print('Interest'.ljust(9),'','End.bal.'.ljust(13))print(''.rjust(6,'-'),'',''.ljust(13,'-'),'',)print(''.rjust(9,'-'),'',''.ljust(9,'-'),'',)print(''.rjust(9,'-'),'',''.ljust(13,'-'),‘‘)defpmt(principal,rate,term):ratePerTwelve=rate/(12*100.0)result=principal*(ratePerTwelve/(1-(1+ratePerTwelve)**(-term)))#Converttodecimalandroundofftotwodecimalplaces.result=Decimal(result)result=round(result,2)returnresult12《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)6.1確定性模型：總回報(bào)plt.figure(figsize=(18,14))#amortization_table(150000,4,180,500)i=0markers=['o','s','D','^','v','*','p','s','D','o','s','D','^','v','*','p','s','D']markersize=[8,8,8,12,8,8,8,12,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8]forextrainrange(100,1700,100):xv,bv,saved=amortization_table(450000,5,360,extra,False)ifextra==0:plt.plot(xv,bv,color=colors[i],lw=2.2,label='Principalonly',marker=markers[i],markersize=markersize[i])else:plt.plot(xv,bv,color=colors[i],lw=2.2,label="Principalplus\$"+str(extra)+str("/month,Saved:\$")+saved,marker=markers[i],markersize=markersize[i])i+=1plt.grid(True)plt.xlabel('Years',fontsize=18)plt.ylabel('MortgageBalance',fontsize=18)plt.title("MortgageLoanFor$350,000WithAdditionalPaymentChart",fontsize=20)plt.legend()plt.show()13《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)6.1確定性模型：總回報(bào)額外儲蓄與本金儲蓄

不同貸款金額的儲蓄

14《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)本金抵押貸款本金的額外金額

6.1確定性模型：總回報(bào)importmatplotlib.pyplotasplt#setthesavingsvaluefrompreviousexampleyvals1=[101000,111000,121000,131000,138000,143000,148000,153000,158000]yvals2=[130000,142000,155000,160000,170000,180000,190000,194000,200000]yvals3=[125000,139000,157000,171000,183000,194000,205000,212000,220000]xvals=[500,600,700,800,900,1000,1100,1200,1300]#initializebubblesthatwillbescaledbubble1=[]bubble2=[]bubble3=[]#scaleitonsomethingthatcanbedisplayed#Itshouldbescaledto1000,butdisplaywillbetoobig#sowechoosetoscaleby5%(dividetheseby20againtorelatetorealvalues)15《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)6.1確定性模型：總回報(bào)foriinrange(0,9):bubble1.append(yvals1[i]/20)bubble2.append(yvals2[i]/20)bubble3.append(yvals3[i]/20)

#plotyvalueswithscaledbybubblesizes.#Ifbubblesarenotscaled,theydon'tfitwell.

fig,ax=plt.subplots(figsize=(10,12))plt1=ax.scatter(xvals,yvals1,c='#d82730',s=bubble1,alpha=0.5)plt2=ax.scatter(xvals,yvals2,c='#2077b4',s=bubble2,alpha=0.5)plt3=ax.scatter(xvals,yvals3,c='#ff8010',s=bubble3,alpha=0.5)

#Setthelabelsandtitleax.set_xlabel('ExtraDollarAmount',fontsize=16)ax.set_ylabel('Savings',fontsize=16)ax.set_title('MortgageSavings(PayingExtraEveryMonth)',fontsize=20)

#setxandylimits

ax.set_xlim(400,1450)

ax.set_ylim(90000,230000)

ax.grid(True)

ax.legend((plt1,plt2,plt3),('$250,000Loan','$350,000Loan','$450,000Loan'),

scatterpoints=1,loc='upperleft',markerscale=0.17,fontsize=10,ncol=1)

fig.tight_layout()

plt.show()16《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)6.1確定性模型：總回報(bào)不同貸款金額額外支付節(jié)省費(fèi)用6.2.1蒙特卡洛模擬6.2.2投資組合模型6.2.3模擬模型6.2.4幾何布朗運(yùn)動模擬6.2.5基于擴(kuò)散的模擬6.2隨機(jī)模型17《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)6.2.1

蒙特卡洛模擬18《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)蒙特卡洛模擬：概率模擬

（1）用于理解風(fēng)險(xiǎn)和不確定性

（2）以隨機(jī)數(shù)作為輸入，迭代進(jìn)行

（3）輸出結(jié)果帶有隨機(jī)性確定性模型示意圖隨機(jī)性模型示意圖

19《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)6.2.1

蒙特卡洛模擬20《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)6.2.1

蒙特卡洛模擬例：利用Python代碼進(jìn)行利潤最大化importnumpyasnpfrommathimportlogimportmatplotlib.pyplotaspltx=[]y=[]#EquationthatdefinesProfitdefgenerateProfit(d):globalsifd>=s:return0.6*selse:return0.6*d-0.4*(s-d)#Althoughycomesfromuniformdistributionin[80,140]#wearerunningsimulationfordin[20,305]maxprofit=0forsinrange(20,305):foriinrange(1,1000):#generatearandomvalueofdd=np.random.randint(10,high=200)profit=generateProfit(d)ifprofit>maxprofit:maxprofit=profitx.append(s)y.append(log(maxprofit))#plottedonlogscaleplt.plot(x,y)print("MaxProfit:",maxprofit)21《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)6.2.1

蒙特卡洛模擬蒙特卡羅模擬利潤值

6.2.1

蒙特卡洛模擬22《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)股票價(jià)格波動率：

（1）利用統(tǒng)計(jì)指標(biāo)（如標(biāo)準(zhǔn)差）衡量波動率

（2）利用金融時(shí)間序列的收盤價(jià)格作為如數(shù)

（3）查看給定時(shí)間范圍內(nèi)的價(jià)格波動率

（4）免費(fèi)平臺Tushare實(shí)時(shí)下載股票交易數(shù)據(jù)例：收益波動率

23《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)6.2.1

蒙特卡洛模擬importtushareastsimportpandasaspd#注冊tushare賬號，可獲得tokents.set_token('yourtoken')pro=_api()

defget_data(code,start,end):df=pro.daily(ts_code=code,start_date=start,end_date=end)df=df.sort_values(by="trade_date")#設(shè)置把日期作為索引df.index=pd.to_datetime(df.trade_date)df=df[['open','high','low','close','vol']]returndfimportnumpyasnpdf=get_data(code="000001.SZ",start="20180101",end="20220101")print(df)r,g,b=(31,119,180)colornow=(r/255.,g/255.,b/255.)print(df)defget_vol(df):df['volatility']=np.log(df['close']/df['close'].shift(1))returndfdf=get_vol(df)df[['close','volatility']].plot(figsize=(12,10),subplots=True,color=colornow)例：收益波動率（多只股票）24《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)6.2.1

蒙特卡洛模擬code_list=["000001.SZ","000021.SZ","000009.SZ","000012.SZ","000004.SZ","000008.SZ"]df_list=[get_vol(get_data(code=code,start="20180101",end="20220101"))forcodeincode_list]foriinrange(len(code_list)):df_list[i]=df_list[i].drop(columns=['open','high','low','close','vol'])df_list[i].columns=[code_list[i]]data=pd.concat(df_list,axis=1,join="inner").dropna()print(data)

importpandasaspdimportmatplotlib.pyplotasplt

#toplotalltogether#vstoxx_short.plot(figsize=(15,14))plt.plot(data,linewidth=0.5,label=data.columns.values)plt.legend()plt.grid()plt.title('SZ’)plt.show()6.2.1

蒙特卡洛模擬25《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)隱含波動率：（1）基于Black-Scholes-Merton模型

（2）估算期權(quán)價(jià)格時(shí)常用（3）隱含波動率：給定歐式看漲期權(quán)C,隱含波動率可由對數(shù)回報(bào)的標(biāo)準(zhǔn)差計(jì)算得到。關(guān)于波動率的期權(quán)定價(jià)公式的偏導(dǎo)數(shù)稱為Vega。該指標(biāo)指示期權(quán)價(jià)格的變動方向以及波動率提高1%的期權(quán)價(jià)格變動的程度

。

例：隱含波動率

26《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)6.2.1

蒙特卡洛模擬frommathimportlog,sqrt,expfromscipyimportstatsimportpandasaspdimportmatplotlib.pyplotasplt

colors=[(31,119,180),(174,199,232),(255,128,0),(255,15,14),(44,160,44),(152,223,138),(214,39,40),(255,152,150),(148,103,189),(197,176,213),(140,86,75),(196,156,148),(227,119,194),(247,182,210),(127,127,127),(199,199,199),(188,189,34),(219,219,141),(23,190,207),(158,218,229)]#ScaletheRGBvaluestothe[0,1]range,whichistheformatmatplotlibaccepts.foriinrange(len(colors)):r,g,b=colors[i]colors[i]=(r/255.,g/255.,b/255.)例：隱含波動率

27《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)6.2.1

蒙特卡洛模擬defblack_scholes_merton(S,r,sigma,X,T):S=float(S)logsoverx=log(S/X)halfsigmasquare=0.5*sigma**2sigmasqrtT=sigma*sqrt(T)d1=(logsoverx+((r+halfsigmasquare)*T))/sigmasqrtTd2=(logsoverx+((r-halfsigmasquare)*T))/sigmasqrtT#stats.norm.cdf—>cumulativedistributionfunctionvalue=S*stats.norm.cdf(d1,0.0,1.0)-X*exp(-r*T)*stats.norm.cdf(d2,0.0,1.0)returnvalue

defvega(S,r,sigma,X,T):S=float(S)logsoverx=log(S/X)halfsigmasquare=0.5*sigma**2sigmasqrtT=sigma*sqrt(T)d1=(logsoverx+((r+halfsigmasquare)*T))/sigmasqrtTvega=S*stats.norm.cdf(d1,0.0,1.0)*sqrt(T)returnvega

defimpliedVolatility(S,r,sigma_est,X,T,Cstar,it):

foriinrange(it):

numer=(black_scholes_merton(S,r,sigma_est,X,T)-Cstar)

denom=vega(S,r,sigma_est,X,T)

sigma_est-=numer/denom

returnsigma_est例：隱含波動率

28《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)6.2.1

蒙特卡洛模擬importpandasaspdfutures_data=pd.read_csv("./data/futures_data.csv")options_data=pd.read_csv("./data/options_data.csv")

options_data['IMP_VOL']=0.0V0=17.6639#theclosingvalueoftheindexr=0.04#riskfreeinterestratesigma_est=2tol=0.5#tolerancelevelformoneynessforoptioninoptions_data.index:#iteratingoveralloptionquotesfutureval=futures_data[futures_data['MATURITY']==options_data.loc[option]['MATURITY']]['PRICE'].values[0]#pickingtherightfuturesvalueif(futureval*(1-tol)<options_data.loc[option]['STRIKE']<futureval*(1+tol)):impliedVol=impliedVolatility(V0,r,sigma_est,options_data.loc[option]['STRIKE'],options_data.loc[option]['TTM'],options_data.loc[option]['PRICE'],#Cnit=100)#iterationsoptions_data['IMP_VOL'].loc[option]=impliedVol

例：隱含波動率

29《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)6.2.1

蒙特卡洛模擬

plot_data=options_data[options_data['IMP_VOL']>0]maturities=sorted(set(options_data['MATURITY']))plt.figure(figsize=(15,10))i=0formaturityinmaturities:data=plot_data[options_data.MATURITY==maturity]#selectdataforthismaturityplot_args={'lw':3,'markersize':9}plt.plot(data['STRIKE'],data['IMP_VOL'],label=maturity,marker='o',color=colors[i],**plot_args)i+=1plt.grid(True)plt.xlabel('Strikerate$X$',fontsize=18)plt.ylabel(r'Impliedvolatilityof$\sigma$',fontsize=18)plt.title('ShortMaturityWindow(VolatilitySmile)',fontsize=22)plt.legend()plt.show()

針對歐洲VSTOXX的罷工率隱含波動率

6.2.2投資組合估值30《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)投資組合估值：

（1）估計(jì)投資組合當(dāng)前價(jià)值

（2）適用于金融資產(chǎn)、金融負(fù)債等

（3）構(gòu)建投資組合、動態(tài)調(diào)整

（4）免費(fèi)平臺Tushare實(shí)時(shí)下載股票交易數(shù)據(jù)例：投資組合估值

31《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)6.2.2投資組合估值importakshareasakimportpandasaspdimportdatetime

defget_fund_data(code,start,end):df=ak.fund_etf_hist_sina(symbol=code)[['date','close']]df=df[df['date']>=datetime.date.fromisoformat(start)]df=df[df['date']<=datetime.date.fromisoformat(end)]df.index=df['date']df=df.drop(columns=['date'])df.columns=[code]returndf

df1=get_fund_data(code="sz169101",start="2018-01-01",end="2022-01-01")df2=get_fund_data(code="sz169102",start="2018-01-01",end="2022-01-01")df3=get_fund_data(code="sz169103",start="2018-01-01",end="2022-01-01")data=pd.concat([df1,df2,df3],axis=1,join="inner")print(data)三只基金收盤價(jià)走勢例：投資組合估值

32《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)6.2.2投資組合估值importnumpyasnp#convertpricestologreturnsretn=data.apply(np.log).diff()#makecorrmatrixretn.corr()#makescatterplottoshowcorrelationpd.plotting.scatter_matrix(retn,figsize=(10,10))plt.show()#somemorestatsretn.skew()retn.kurt()三只基金收益相關(guān)性圖

33《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)6.2.3模擬模型模擬仿真模型：

（1）基于隨機(jī)數(shù)和概率

（2）對真實(shí)數(shù)據(jù)和世界進(jìn)行仿真實(shí)驗(yàn)

（3）科研和實(shí)踐的試錯和檢驗(yàn)工具

（4）復(fù)雜計(jì)算中的重要估算工具 ……幾何布朗運(yùn)動定義：

34《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)6.2.4幾何布朗運(yùn)動模擬

例：幾何布朗運(yùn)動模擬

35《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)importmatplotlib.pyplotaspltimportnumpyasnprect=[0.1,5.0,0.1,0.1]fig=plt.figure(figsize=(10,10))T=2mu=0.1sigma=0.04S0=20dt=0.01N=round(T/dt)t=np.linspace(0,T,N)#StandarenormaldistribW=np.random.standard_normal(size=N)W=np.cumsum(W)*np.sqrt(dt)X=(mu-0.5*sigma**2)*t+sigma*W#BrownianMotionS=S0*np.exp(X)plt.plot(t,S,lw=2)plt.xlabel("Timet",fontsize=16)plt.ylabel("S",fontsize=16)plt.title("GeometricBrownianMotion(Simulation)",fontsize=18)plt.show()幾何布朗運(yùn)動模擬圖

6.2.4幾何布朗運(yùn)動模擬例：幾何布朗運(yùn)動模擬股票價(jià)格

36《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)6.2.4幾何布朗運(yùn)動模擬importpylab,randomclassStock(object):def__init__(self,price,distribution):self.price=priceself.history=[price]self.distribution=distributionself.lastChange=0defsetPrice(self,price):self.price=priceself.history.append(price)defgetPrice(self):returnself.pricedefwalkIt(self,marketBias,mo):oldPrice=self.pricebaseMove=self.distribution()+marketBiasself.price=self.price*(1.0+baseMove)ifmo:self.price=self.price+random.gauss(.5,.5)*self.lastChangeifself.price<0.01:self.price=0.0self.history.append(self.price)self.lastChange=oldPrice-self.pricedefplotIt(self,figNum):pylab.figure(figNum)pylab.plot(self.history)pylab.title('ClosingPriceSimulationRun-'+str(figNum))pylab.xlabel('Day')pylab.ylabel('Price')例：幾何布朗運(yùn)動模擬股票價(jià)格

37《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)6.2.4幾何布朗運(yùn)動模擬deftestStockSimulation():defrunSimulation(stocks,fig,mo):mean=0.0forsinstocks:fordinrange(numDays):s.walkIt(bias,mo)s.plotIt(fig)mean+=s.getPrice()mean=mean/float(numStocks)pylab.axhline(mean)pylab.figure(figsize=(12,12))numStocks=20numDays=400stocks=[]bias=0.0mo=Falsestartvalues=[100,500,200,300,100,100,100,200,200, 300,300,400,500,00,300,100,100,100,200,200,300]foriinrange(numStocks):volatility=random.uniform(0,0.2)d1=lambda:random.uniform(-volatility,volatility)stocks.append(Stock(startvalues[i],d1))runSimulation(stocks,1,mo)

testStockSimulation()pylab.show()幾何布朗運(yùn)動模擬圖

平方根擴(kuò)散模型，用于對平均回歸量建模（如利率和波動性）

38《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)6.2.5基于擴(kuò)散的模擬

代碼實(shí)例

39《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)6.2.5基于擴(kuò)散的模擬importnumpyasnpimportmatplotlib.pyplotaspltimportnumpy.randomasnprS0=100#initialvaluer=0.05sigma=0.25T=2.0x0=0k=1.8theta=0.24i=100000M=50dt=T/Mdefsrd_euler():xh=np.zeros((M+1,i))x1=np.zeros_like(xh)xh[0]=x0x1[0]=x0fortinrange(1,M+1):xh[t]=(xh[t-1]+k*(theta-np.maximum(xh[t-1],0))*dt+sigma*np.sqrt(np.maximum(xh[t-1],0))*np.sqrt(dt)*npr.standard_normal(i))x1=np.maximum(xh,0)returnx1x1=srd_euler()代碼實(shí)例

40《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)6.2.5基于擴(kuò)散的模擬平方根擴(kuò)散模型plt.figure(figsize=(10,6))plt.hist(x1[-1],bins=30,color='#98DE2f',alpha=0.85)plt.xlabel('value')plt.ylabel('frequency')plt.grid(False)plt.figure(figsize=(12,10))plt.plot(x1[:,:10],lw=2.2)plt.title("Square-RootDiffusion-Simulation")plt.xlabel('Time',fontsize=16)plt.ylabel('IndexLevel',fontsize=16)#plt.grid(True)plt.show()6.3閾值模型-謝林隔離模型41《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)謝林隔離模型（Schelling’s隔離模型）

（1）應(yīng)用場景一:通過在棋盤上放置硬幣并根據(jù)各種規(guī)則移動它們來進(jìn)行實(shí)驗(yàn)。在該實(shí)驗(yàn)中，一個棋盤被用來類比城市，一個棋盤正方形對應(yīng)一個住所，一個正方形對應(yīng)一個社區(qū)。pennies和dimes（視覺上不同）可以代表兩組的吸煙者、非吸煙者、男性、女性、高管、非高管、學(xué)生或教師。模擬規(guī)則指定當(dāng)沒有人從當(dāng)前位置移動則終止，即代表所有人都滿意。若有人不高興，他們就會搬家，即產(chǎn)生隔離。

（2）應(yīng)用場景二：模擬教室選擇的過程，該模型表明即使對相鄰?fù)瑢W(xué)的偏好較弱，也可能出現(xiàn)隔離模式。

（3）該模型類似于動態(tài)聚類過程。6.3閾值模型-謝林隔離模型42《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)謝林隔離模型舉例：

假設(shè)一所高中每種類型有250名學(xué)生，并根據(jù)他們的第一優(yōu)先級將學(xué)生類別分為三種類型：運(yùn)動、高級水平學(xué)術(shù)和常規(guī)，每種類型分別為0、1和2來代表。

這些學(xué)生都住在一個單位廣場上（這里可以想象成一棟教學(xué)樓建筑）。它的位置只是一個點(diǎn)(x,y)，其中0<x,y<1。如果它的12個最近鄰居中有一半或更多屬于同一類型（根據(jù)歐幾里得距離），則它很高興。

每個同學(xué)初始位置是從一個二元均勻分布中獨(dú)立抽取的。我們考慮最終形成的自然隔離結(jié)果。

6.3閾值模型-謝林隔離模型43《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)fromrandomimportuniform,seedfrommathimportsqrtimportmatplotlib.pyplotaspltnum=250#ThesemanyagentsofaparticulartypenumNeighbors=12#NumberofagentsregardedasneighborsrequireSameType=8#Atleastthismanyneighborstobesametypeseed(10)#forreproduciblerandomnumbers6.3閾值模型-謝林隔離模型44《數(shù)據(jù)分析與可視化》第五章統(tǒng)計(jì)學(xué)與機(jī)器學(xué)習(xí)classStudentAgent:def__init__(self,type):#Studentsofdifferenttypewillbeshownincolorsself.type=typeself.show_position()

defshow_position(self):#positionchangedbyusinguniform(x,y)self.position=uniform(0,1),uniform(0,1)

defget_distance(self,other):#returnseuclideandistance