1、模擬程序案例例1,在做拋擲一枚質地均勻的硬幣的試驗中發“正面朝上”的事件(用 1表示)和“正面朝下”的事件A (用0表示)的情況。歷史上一些學者得到的具體試驗結果如下：現在需要利用eviews來模擬上述三位學者的實驗。算法分析：上述三學者的實驗均為二項分布的實驗，可以直接利用eviews產生二項分布隨機數的函數rbinom (n, p).編程如下：workfile binom u 1 2048 series result for !i=0 to 1 smpl 1 2048 series x x(1)=0for !cou=1 to 2048 x(!cou)=rbinom(!i,0.5) next

2、 nextSeries: XSample 1 2048Observations 2048Mean0.507813Median1.000000Maximum1.000000Minimum0.000000Std. Dev.0.500061Skewness-0.031254Kurtosis1.000977Jarque-Bera341.3334Probability0.000000 x.hist 1,2001,000 -800 _600 -400 -200 -0-,1,11,10.00.10.20.30.40.50.60.70.80.91.0例2 (投擲骰子)(1)投擲一顆質地均勻的骰子，令 X表示其

3、出現的點數，分析各點數出 現的頻率的穩定性及變化規律；(2)利用統計的方法，根據“頻率的穩定性”規律求投擲一 枚質地不均勻的骰子出現某點數的概率；(3)演示隨機變量X的數學期望的統計意義。算法分析：根據逆變換法產生來自分布函數 F (x)的隨機數，就要求出F-1 (y),其中F-1 (y) =infx： F (x)y.0&y& 1.質地均勻的骰子各點數出現的頻率的分布函數是F (x) =p (xx) = (i=1) /6, i-1 xy.0y 1=i-1 , (i-1) /6yi/6, i=1 ,，6因而，可先由產生均勻分布隨機數的函數runif (0, 1)抽取y值，再來計算F-1 (y)值

4、即可。程序實現：workfile binom u 1 1000smpl 1 1000series xseries yseries a1series a2series a3series a4series a5series a6for !i=1 to 1000a1(!i)=1/6a2(!i)=2/6a3(!i)=3/6a4(!i)=4/6a5(!i)=5/6a6(!i)=1x(!i)=runif(0,1)if x(!i)=a1(!i) and x(!i)=a2(!i) and x(!i)=a3(!i) and x(!i)=a4(!i) and x(!i)=a5(!i) and x(!i)a6(!i

5、) then y(!i)=6 else y(!i)=7endifendifendifendifendifendifnexty.histSeries: YSample 1 1000Observations 1000Mean Median Maximum Minimum Std. Dev. Skewness KurtosisJarque-Bera Probability3.5240003.0000006.0000001.0000001.680550 -0.002354 1.76118863.944890.0000001.通過已知總體模型得到多組樣本數據，進行多次回歸，驗證回歸結果的特征、性質 最小

6、二乘法的無偏性workfile mc u 1 10vector(10) v1vl.fill 80, 100,120,140,160,180,200,220,240,260mtos(v1,x)!b1=25!b2=0.5matrix(100,2) ffor !k=1 to 100series u=3*nrndseries y=!b1+!b2*x+uequation eq.ls y=c(1)+c(2)* xf(!k,1)=c(1)f(!k,2)=c(2)nextshow fexpand 1 100smpl 1 100mtos(f,gr)freeze ser01.qqplotfreeze ser01.

7、histfreeze ser02.qqplotfreeze ser02.histmatrix(1,2) mm(1,1)=mean(ser01)m(1,2)=mean(ser02)show ma ELONfo sownnauQ.46.46.47 .48.49.50.51.52.53.54.55a ELONf o54321098755555544432161618202224262830308 6 42 2 222o 82 1Quantiles of SER02Quantiles of SER01蒙特卡洛模擬程序：(最終調試成功)store monte carle results in a ser

8、ieschecked 4/1/2004set workfile range to number of monte carle replicationswfcreate mcarle u 1 100create data series for xnote: x is fixed in repeated samplesonly first 10 observations are used (remaining 90 obs missing) series xx.fill 80,100,120,140,160,180,200,220,240,260set true parameter values!

9、beta1=2.5!beta2=0.5set seed for random number generatorrndseed 123456assign number of replications to a control variable!reps=100begin loopfor !I=1 to !repsset sample to estimation samplesmpl 1 10simulate y data (only for 10 obs)series y=!beta1 + !beta2*x+3*nrndregress y on a constant and xEquation

10、eq1.ls y c xset sample to one observationsmpl !I !iand store each coefficient estimate in a seriesseries b1=eq1.coefs(1) series b2=eq1.coefs(2) next end of loopset sample to full samplesmpl 1 100show kernel density eatimate for each coef freeze(gra1) b1.distplot kernel!beta1!beta2drow vertical dashl

11、ine at true parameter value gra1.draw(dashline,bottom,rgb(156,156,156) show gra1freeze(gra2) b2.distplot kerneldraw vertical dashline at true parameter value gra2.draw(dashline,bottom,rgb(156,156,156) show gra2一元回歸參數的分布:Subroutine moni (scalar n, scalar sum, scalar param1, scalar param2, scalar type

12、) For !m=1 to nX(!m)=rnd*100NextFor !n=1 to sumIf type=0 thenFor !m=1 to nU(!m)=nrndNextEndifIf type=1 thenFor !m=1 to nIf rnd0.5 thenU(!m)=rndElseU(!m)=rnd*(-1)EndifNextEndifIf type=2 thenU(1)=nrndFor !m=2 to nU(!m)=u(!m-1)+nrndNextEndifIf type=3 then For !m=1 to nIf rnd=2freeze(statby_%1) %1.statb

13、y(nomean,nostd)cao_%1nextfor %1 r0 b0 t0 dw0 r20 f0 r1 b1 t1 dw1 r21 f1 r2 b2 t2 dw2 r22 f2freeze(hist_%1) %1.histnext偽回歸相關系數模擬(不行)workfile corr u 1 500 series resultfor !i=1 to 500smpl 1 100series x=nrnd series y=nrnd series xx series yyscalar sum1=0scalar sum2=0 for !counter=1 to 100 sum1=sum1+x(!

14、counter) sum2=sum2+y(!counter) xx(!counter)=sum1 yy(!counter)=sum2 next nextscalar r=cor(xx,yy) result(!i)=r時變貝塔系數的模擬(執行不了)：BETA.PRG (3/7/2007) Time varying beta demonstrates several ways to obtain beta between assets1) constant beta2) rolling beta by regression/moving cov/var3) using state space4)

15、using multivariate ARCHChecked 3/20/2007change path to program path%path = runpathcd %pathload workfileload 僅.wf1dependent variables of series must be continuous smpl allseries y1 = pch(index)series y2 = pch(jy) calculate the constant beta using OLSsmpl 1990 lastequation constant_beta.ls y2 c y1seri

16、es beta_const=c(2)calculating time varying beta with rolling regression for a bi-variate case can use moving cov/var instead of OLS regression!ssize = 200series beta_roll=movcov(y1,y2,!ssize)/movvar(y1,!ssize)code for running a rolling regression:commented out right now!length = obs(y1)equation roll

17、_beta.ls y2 c y1show roll_betafor !i = 1 to !length-!ssize+1smpl first+!i-1 first+!i+!ssize-2equation roll_beta.ls y2 c y1smpl first+!i+!ssize-2 first+!i+!ssize-2beta_roll = roll_beta.coefs(2)nextcalculate beta with State Spacevia a time-varying coefficient for Y1 smpl 1990 lastsspace ssbetassbeta.a

18、ppend y2=c(1)+sv1*y1+var=exp(c(2)ssbeta.append state sv1 = sv1(-1)ssbeta.mlssbeta.makestates beta_* rename beta_sv1 beta_sscalculate beta with system ARCHby estimating the covariance and variance of the two series using Multivariate ARCHsystem arbetaarbeta.append y1 = c(1)arbeta.append y2 = c(2)arbeta.arch Diagvech c(indef) arch(1,indef) garch(1,indef)arbeta.makegarch(name=arch)series beta_arch = arch01_02/arch01 display the different betasgroup betas_ls_roll beta_const beta_rollgroup betas_roll_ss_arch beta_roll beta_ss