蒙特卡罗 模拟实验

蒙特卡罗模拟循环;
怀特检验循环100次的code
x=c(4.36,4.50,4.50,4.51,4.55,4.62,4.68,4.73,4.75,4.8,4.92,4.95,5.1,5.23,5.36,5.38,5.39,5.51,5.65,5 .84,5.86,5.92,6.23,6.59,6.91,6.96,7.67,8.48,9.21,9.24,10.99)
f=function(){
a=rnorm(31,0,log(x))
y=0.8*x+a-0.2
res1=lm(y~x)
ei_square=residuals(res1)^2
res2=lm(ei_square~x+I(x^2))
ts=summary(res2)$"r.squared"*length(x)
}
res=c()
for(i in 1:1000) res[i]=f()
qchisq(0.95,2)
5.991465
pvalue=mean(res>= 5.991465)
pvalue
park 检验蒙特卡罗模拟循环
x=c(4.36,4.50,4.50,4.51,4.55,4.62,4.68,4.73,4.75,4.8,4.92,4.95,5.1,5.23,5.36,5.38,5.39,5.51,5.65,5 .84,5.86,5.92,6.23,6.59,6.91,6.96,7.67,8.48,9.21,9.24,10.99)
f=function(){
a=rnorm(31,0, log(x))
y=0.8*x+a-0.2
res1=lm(y~x)
ei_square=residuals(res1)^2
res2=lm(log(ei_square)~log(x))
fs=summary(res2)$"fstatistic"
}
res=c()
for(i in 1:1000) res[i]=f()
res
qf(0.95,1,29)
pvalue=mean(res>=4.182964)
pvalue
###############################################
glzezser 检验
x=c(4.36,4.50,4.50,4.51,4.55,4.62,4.68,4.73,4.75,4.8,4.92,4.95,5.1,5.23,5.36,5.38,5.39,5.51,5.65,5 .84,5.86,5.92,6.23,6.59,6.91,6.96,7.67,8.48,9.21,9.24,10.99)
f=function(){
a=rnorm(31,0,x)
y=0.8*x+a-0.2
res1=lm(y~x)
ei=residuals(res1)
res_GL=lm(abs(ei)~I(x^-1)) 当指数为其他时，类似只需改变I(x^h)里面的h值即可) fs=summary(res_GL)$"f..statistic" #得出回归方程，检验回归系数，分析比较结果。

}
res=c()
for(i in 1:100) res[i]=f()
res
qf(0.95,1,29)
qf
pvalue=mean(res>=4.182964)
pvalue
G-Q检验
x=c(4.36,4.50,4.50,4.51,4.55,4.62,4.68,4.73,4.75,4.8,4.92,4.95,5.1,5.23,5.36,5.38,5.39,5.51,5.65,5 .84,5.86,5.92,6.23,6.59,6.91,6.96,7.67,8.48,9.21,9.24,10.99)
f=function(){
a=rnorm(31,0,log(x))
y=0.8*x+a-0.2
res1=lm(y~x)
ei=residuals(res1)
x.sorted=sort(x)
resG1=lm(y[1:11]~x.sorted[1:11])
resG2=lm(y[21:31]~x.sorted[21:31])
GQ_TEST=sum(residuals(resG1)^2)/sum(residuals(resG2)^2)
fs=summary(GQ_TEST)$"fstatistic"
}
res=c()
for(i in 1:1000) res[i]=f()
res
qf(0.95,9,9)
qf
pvalue=mean(res>=3.178893)
pvalue
Breusch-Pagan/Godfre检验
x=c(4.36,4.50,4.50,4.51,4.55,4.62,4.68,4.73,4.75,4.8,4.92,4.95,5.1,5.23,5.36,5.38,5.39,5.51,5.65,5 .84,5.86,5.92,6.23,6.59,6.91,6.96,7.67,8.48,9.21,9.24,10.99)
f=function(){
a=rnorm(31,0,log(x))
y=0.8*x+a-0.2
res1=lm(y~x)
ei=residuals(res1)
T=30
reps=numeric(T)
agm=sum(ei_square[1])/length(x)
ei_square=residuals(res1)^2
for(i in 2:31){reps[i]=sum(ei_square[c(1:i)])/length(x)}
zgmi=c(agm,reps)[-2]
p=ei_square/zgmi
fs=summary(p)$"fstatistic"
}
res=c()
for(i in 1:1000) res[i]=f()
res
qf(0.95,1,29)
qf
pvalue=mean(res>=4.182964)
pvalue
park
x=c(4.36,4.50,4.50,4.51,4.55,4.62,4.68,4.73,4.75,4.8,4.92,4.95,5.1,5.23,5.36,5.38,5.39,5.51,5.65,5 .84,5.86,5.92,6.23,6.59,6.91,6.96,7.67,8.48,9.21,9.24,10.99)
f=function(){
a=rnorm(31,0,x^2)
y=0.8*x+a-0.2
res1=lm(y~x)
ei_square=residuals(res1)^2
res2=lm(log(ei_square)~log(x))
fs=summary(res2)$"fstatistic"
}
res=c()
for(i in 1:1000) res[i]=f()
res
qf(0.95,1,29)
pvalue=mean(res>=4.182964)
pvalue
glzezser
x=c(4.36,4.50,4.50,4.51,4.55,4.62,4.68,4.73,4.75,4.8,4.92,4.95,5.1,5.23,5.36,5.38,5.39,5.51,5.65,5 .84,5.86,5.92,6.23,6.59,6.91,6.96,7.67,8.48,9.21,9.24,10.99)
f=function(){
a=rnorm(31,0,x)
y=0.8*x+a-0.2
res1=lm(y~x)
ei=residuals(res1)
res_GL=lm(abs(ei)~I(x^-1))
fs=summary(res_GL)$"fstatistic"
}
res=c()
for(i in 1:100) res[i]=f()
res
qf(0.95,1,29)
qf
pvalue=mean(res>=4.182964)
Pvalue。