模拟程序案例

例 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)=0

for !cou=1 to 2048

x(!cou)=@rbinom(!i,0.5)

next

next

x.hist

word/media/image1_1.png

例 2（投掷骰子）（1）投掷一颗质地均匀的骰子，令 X 表示其出现的点数，分析各点数出现的频率的稳定性及变化规律；（2）利用统计的方法，根据“频率的稳定性”规律求投掷一枚质地不均匀的骰子出现某点数的概率；（3）演示随机变量 X 的数学期望的统计意义。

算法分析：根据逆变换法产生来自分布函数F（x）的随机数，就要求出 F-1（y），其中 F-1（y）＝inf{x：F（x）≥y}.0≤y≤1.质地均匀的骰子各点数出现的频率的分布函数是

F（x）=p（x≥x）=（i=1）/6，i-1≤x＜i，i=1，2，…，7

可求得

F-1（y）=inf{x:F（x）≥y}.0≤y≤1=i-1，（i-1）/6≤y＜i/6，i=1，…，6

因而，可先由产生均匀分布随机数的函数 @runif（0，1）抽取 y 值，再来计算 F-1（y）值即可。

程序实现：

workfile binom u 1 1000

smpl 1 1000

series x

series y

series a1

series a2

series a3

series a4

series a5

series a6

for !i=1 to 1000

a1(!i)=1/6

a2(!i)=2/6

a3(!i)=3/6

a4(!i)=4/6

a5(!i)=5/6

a6(!i)=1

x(!i)=@runif(0,1)

if x(!i)

else if x(!i)>=a1(!i) and x(!i)

else if x(!i)>=a2(!i) and x(!i)(!i) then y(!i)=3

else if x(!i)>=a3(!i) and x(!i)(!i) then y(!i)=4

else if x(!i)>=a4(!i) and x(!i)(!i) then y(!i)=5

else if x(!i)>=a5(!i) and x(!i)(!i) then y(!i)=6

else y(!i)=7

endif

endif

endif

endif

endif

endif

next

y.hist

word/media/image2_1.png

1.通过已知总体模型得到多组样本数据，进行多次回归，验证回归结果的特征、性质

最小二乘法的无偏性

workfile mc u 1 10

vector(10) v1

v1.fill 80, 100,120,140,160,180,200,220,240,260

mtos(v1,x)

!b1=25

!b2=0.5

matrix(100,2) f

for !k=1 to 100

series u=3*nrnd

series y=!b1+!b2*x+u

equation eq.ls y=c(1)+c(2)* x

f(!k,1)=c(1)

f(!k,2)=c(2)

next

show f

expand 1 100

smpl 1 100

mtos(f,gr)

freeze ser01.qqplot

freeze ser01.hist

freeze ser02.qqplot

freeze ser02.hist

matrix(1,2) m

m(1,1)=@mean(ser01)

m(1,2)=@mean(ser02)

show m

word/media/image3_1.pngword/media/image4_1.png

蒙特卡洛模拟程序：(最终调试成功)

'store monte carle results in a series

'checked 4/1/2004

'set workfile range to number of monte carle replications

wfcreate mcarle u 1 100

'create data series for x

'note: x is fixed in repeated samples

'only first 10 observations are used (remaining 90 obs missing)

series x

x.fill 80,100,120,140,160,180,200,220,240,260

'set true parameter values

!beta1=2.5

!beta2=0.5

'set seed for random number generator

rndseed 123456

'assign number of replications to a control variable

!reps=100

'begin loop

for !I=1 to !reps

'set sample to estimation sample

smpl 1 10

'simulate y data (only for 10 obs)

series y=!beta1+!beta2*x+3*nrnd

'regress y on a constant and x

Equation eq1.ls y c x

'set sample to one observation

smpl !I !i

'and store each coefficient estimate in a series

series b1=eq1.@coefs(1)

series b2=eq1.@coefs(2)

next

'end of loop

'set sample to full sample

smpl 1 100

'show kernel density eatimate for each coef

freeze(gra1) b1.distplot kernel

'drow vertical dashline at true parameter value

gra1.draw(dashline,bottom,rgb(156,156,156)) !beta1

show gra1

freeze(gra2) b2.distplot kernel

'draw vertical dashline at true parameter value

gra2.draw(dashline,bottom,rgb(156,156,156)) !beta2

show gra2

一元回归参数的分布：

Subroutine moni (scalar n, scalar sum, scalar param1, scalar param2, scalar type)

For !m=1 to n

X(!m)=@rnd*100

Next

For !n=1 to sum

If type=0 then

For !m=1 to n

U(!m)=@nrnd

Next

Endif

If type=1 then

For !m=1 to n

If @rnd<0.5 then

U(!m)=@rnd

Else

U(!m)=@rnd*(-1)

Endif

Next

Endif

If type=2 then

U(1)=@nrnd

For !m=2 to n

U(!m)=u(!m-1)+@nrnd

Next

Endif

If type=3 then

For !m=1 to n

If @rnd<0.5 then

U(!m)=@nrnd

Else

U(!m)=@nrnd*2

Endif

Next

Endif

For !m=2 to n

Y(!m)=param1+param2*x(!m)+u(!m)

Next

Equation eq.ls y c x

B1(!n)=eq.@coefs(1)

B2(!n)=eq.@coefs(2)

Next

Endsub

Workfile moni u 1 10000

Series b1

Series b2

Series x

Series y

Series u

Call moni(100,10000,20,0.8,3)

B1.hist

B2.hist

word/media/image5_1.pngword/media/image6_1.png

生成季度虚拟变量的程序：

wfcreate dumtest q 1970 1990

%start = "1972:1"

%end = "1979:4"

for !i = @dtoo(%start) to @dtoo(%end)

%obsstr = @otod(!i)

if (@mid(%obsstr, 5, 1) = ":") then

%name = "d_" + @left(%obsstr, 4) + "_" + @mid(%obsstr, 6)

else

%name = "d_" + %obsstr

endif

smpl @all

series {%name} = 0

smpl {%obsstr} {%obsstr}

series {%name} = 1

next

下面这个执行不了：

!n=10000

wfcreate cao u 1 !n

matrix(!n,18) m

group g

for %1 r0 b0 t0 dw0 r20 f0 r1 b1 t1 dw1 r21 f1 r2 b2 t2 dw2 r22 f2

series {%1}=0

g.add {%1}

next

for !k=1 to !n

smpl 1 122

for %1 y0 x0 y1 x1 y2 x2

series {%1}=@rnorm

next

smpl 2 122

x1=x1(-1)+@rnorm

y1=y1(-1)+@rnorm

smpl 2 122

x2=2*x2(-1)-x2(-2)+@rnorm

y2=2*y2(-1)-y2(-2)+@rnorm

smpl 1 122

for %1 %2 %3 1 y0 x0 7 y1 x1 13 y2 x2

m(!k,{%1})=@cor({%2},{%3})

equation eq.ls {%2} c {%3}

m(!k,{%1}+1)=eq.c(2)

m(!k,{%1}+2)=eq.@tstats(2)

m(!k,{%1}+3)=eq.@dw

m(!k,{%1}+4)=eq.@r2

m(!k,{%1}+5)=eq.@f

next

next

smpl 1 !T

mtos(m,g)

for %1 t0 t1 t2

series cao_{%1}=@abs({%1})>=2

freeze(statby_{%1}) {%1}.statby(nomean,nostd)cao_{%1}

next

for %1 r0 b0 t0 dw0 r20 f0 r1 b1 t1 dw1 r21 f1 r2 b2 t2 dw2 r22 f2

freeze(hist_{%1}) {%1}.hist

next

伪回归相关系数模拟（不行）

workfile corr u 1 500

series result

for !i=1 to 500

smpl 1 100

series x=nrnd

series y=nrnd

series xx

series yy

scalar sum1=0

scalar sum2=0

for !counter=1 to 100

sum1=sum1+x(!counter)

sum2=sum2+y(!counter)

xx(!counter)=sum1

yy(!counter)=sum2

next

next

scalar r=@cor(xx,yy)

result(!i)=r

时变贝塔系数的模拟（执行不了）：

' BETA.PRG (3/7/2007)

' Time varying beta

' demonstrates several ways

' to obtain beta between assets

' 1) constant beta

' 2) rolling beta by regression/moving cov/var

' 3) using state space

' 4) using multivariate ARCH

' Checked 3/20/2007

'change path to program path

%path = @runpath

cd %path

' load workfile

load fx.wf1

' dependent variables of series must be continuous

smpl @all

series y1 = @pch(index)

series y2 = @pch(jy)

'------------------------------------------------------------------------

' calculate the constant beta using OLS

'------------------------------------------------------------------------

smpl 1990 @last

equation constant_beta.ls y2 c y1

series 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 = 200

series beta_roll=@movcov(y1,y2,!ssize)/@movvar(y1,!ssize)

' code for running a rolling regression:

' commented out right now

'!length = @obs(y1)

'equation roll_beta.ls y2 c y1

'show roll_beta

'for !i = 1 to !length-!ssize+1

' smpl @first+!i-1 @first+!i+!ssize-2

' equation roll_beta.ls y2 c y1

' smpl @first+!i+!ssize-2 @first+!i+!ssize-2

' beta_roll = roll_beta.@coefs(2)

'next

'------------------------------------------------------------------------

' calculate beta with State Space

' via a time-varying coefficient for Y1

'------------------------------------------------------------------------

smpl 1990 @last

sspace ssbeta

ssbeta.append y2=c(1)+sv1*y1+[var=exp(c(2))]

ssbeta.append @state sv1 = sv1(-1)

ssbeta.ml

ssbeta.makestates beta_*

rename beta_sv1 beta_ss

'------------------------------------------------------------------------

' calculate beta with system ARCH

' by estimating the covariance and variance of

' the two series using Multivariate ARCH

'------------------------------------------------------------------------

system arbeta

arbeta.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 betas

'------------------------------------------------------------------------

group betas_ls_roll beta_const beta_roll

group betas_roll_ss_arch beta_roll beta_ss beta_arch

show betas_ls_roll.line

show betas_roll_ss_arch.line

在生成服从某种分布（标准正态分布白噪声序列）的随机序列的基础上 生成各种ARIMA序列：AR(1)、MA(1)、ARMA(1,1)、ARIMA(1,1,1)

⏹ Workfile random2 u 1 1000

⏹ Series u=@nrnd

⏹ U(1)=0

⏹ Series x1

⏹ X1(1)=0

⏹ Series x2

⏹ X2(1)=0

⏹ Series x3

⏹ X3(1)=0

⏹ Series x4

⏹ X4(1)=0

⏹ Smpl 2 1000

⏹ x1=0.8*x1(-1)+u

⏹ X2=u+0.8*u(-1)

⏹ X3=0.8*x3(-1)+u+0.8*u(-1)

⏹ X4=x4(-1)+x3

⏹ Smpl @all

⏹ group my x1 x2 x3

⏹ show my.line

word/media/image7_1.png

⏹ workfile simu u 1 1000

⏹ series x1=@runif(0,1)

⏹ series x2=@runif(0,1)

⏹ series x3=@runif(0,1)

⏹ group g x1 x2 x3

⏹ g.line

⏹ g.bar

⏹ scat x1 x2

⏹ sort x1

⏹ genr y=10+2*x1

⏹ y.line

⏹ freeze y

⏹ y.hist

⏹ scalar yy=@mean(y)

⏹ scalar yy1=@skew(y)

⏹ scalar yy3=@kurt(y)

word/media/image8_1.png

eviews中的蒙特卡洛模拟程序