问题是这样,这是一个7元7次方程组,有四个外生变量是我关注的,然后运行其中的每一个变量,都会导致内生变量的解有单调性。现在想同时看4个参数分别从0.1取到0.8,等于说最后应该有8x8x8x8个结果,然后我在其中选个最大的,可是现在程序运行不出来,好像因为初始值的取值有问题,请各位高手帮帮忙,多谢了!!! 我用的是[R]程序。
rm(list=ls(all=TRUE))
library(numDeriv)
g=function(sigma,del,baby,sister)
{
init=c(0.9,0.9,0.7,0.1,0.1,1,1)
pbar=1
beta=0.72
r=0.0033
alpha=0.28
v=1
lambda=0.042
m=0.51
c=0.42
eta=0.4
t=3
delta=1
ito=0.95
f=0
par=optim(init,function(x)
{
y=c(0.9,0.9,0.7,0.1,0.1,1,1)
y[1]=v*x[1]-baby-x[4]+eta*(x[7]^t)/t-(x[7]*m*(x[3]^alpha)*beta*v*(pbar-x[1]))/(r+lambda)+(del*m*(x[3]^(alpha))*beta*v*(x[7]*(1-x[1])-x[6]+x[6]*x[2]))/((r+del)*(r+lambda))+((del*(baby-sister-(eta*(x[7]^t-x[6]^t)/t)))/((r+del)))+(lambda*v/((r+lambda)*pbar))*((pbar^2-x[1]^2)/2-(pbar-x[1])*x[1])
y[2]=v*x[2]-sister-x[5]+(eta*(x[6]^t)/t)-(x[6]*m*(x[3]^(alpha))*beta*v*(pbar-x[2]))/(r+lambda)+(sigma*m*(x[3]^(alpha))*beta*v*(x[7]*(1-x[1])-x[6]+x[6]*x[2]))/((r+del)*(r+lambda))+sigma*(baby-sister-(eta*(x[7]^t-x[6]^t)/t))/((r+del))+(beta*sigma*v*(x[2]-x[1]))/(r+lambda)+(lambda*v/((r+lambda)*pbar))*((pbar^2-x[2]^2)/2-(pbar-x[2])*x[2])
y[3]=c/((1-beta))-m*(x[3]^(alpha-1))*((((x[3]^(alpha))*x[6]*m*v*(pbar-x[1])/(((lambda*x[2]*del/sigma)+del+x[6]*m*(x[3]^(alpha)))*(r+lambda))+(lambda*x[2]*del/sigma+del)*v*(pbar-x[2])/((lambda*x[2]*del/sigma+del+x[6]*m*(x[3]^(alpha)))*(r+lambda)))))
y[4]=x[4]*((del+x[7]*m*(x[3]^(alpha)))/(lambda*x[1]))*(1/(del/sigma+(lambda*x[2]*del/sigma+del)/(x[6]*m*(x[3]^(alpha)))+(del+x[7]*m*(x[3]^(alpha)))/(lambda*x[1])+1))+x[5]*(del/sigma)*(1/(del/sigma+(lambda*x[2]*del/sigma+del)/(x[6]*m*(x[3]^(alpha)))+(del+x[7]*m*(x[3]^(alpha)))/(lambda*x[1])+1))-baby*(1/(del/sigma+(lambda*x[2]*del/sigma+del)/(x[6]*m*(x[3]^(alpha)))+(del+x[7]*m*(x[3]^(alpha)))/(lambda*x[1])+1))-sister*((lambda*x[2]*del/sigma+del)/(x[6]*m*(x[3]^(alpha))))*(1/(del/sigma+(lambda*x[2]*del/sigma+del)/(x[6]*m*(x[3]^(alpha)))+(del+x[7]*m*(x[3]^(alpha)))/(lambda*x[1])+1))
y[5]=x[5]*((beta*v*0.5*(1-x[1])+baby+x[4]-eta*(x[7]^t)/t+(beta*v*(1-x[1])*x[7]*m*(x[3]^(alpha))/(r+lambda))-(del*beta*v*(x[7]-x[7]*x[1]-x[6]+x[6]*x[2])*m*(x[3]^(alpha))/((r+del)*(r+lambda)))-(lambda*beta*v*((1-x[1])^2)*0.5/(r+lambda))-(del*(baby-sister-(eta*(x[7]^t-x[6]^t)/t))/(r+del))))-x[4]*((beta*v*0.5*(1-x[2])-eta*(x[6]^t)/t-(lambda*beta*v*((1-x[2])^2)^0.5/(r+lambda))-(sigma*beta*v*(x[7]-x[7]*x[1]-x[6]+x[6]*x[2])*((m*(x[3]^(alpha))/(r+del))))/(r+lambda)-sigma*beta*(x[2]-x[1])/(r+lambda)+sister+x[5]-(sigma*(baby-sister-(eta*(x[7]^t-x[6]^t)/t))/(r+del))+beta*x[6]*m*(x[3]^(alpha))*v*(1-x[2])/(r+lambda)))
y[6]=x[6]-(m*(x[3]^(alpha))*beta*(1-x[2])/(eta*(r+lambda)))^(1/(t-1))
y[7]=x[7]-x[6]*(((1-x[1])/(1-x[2]))^(1/(t-1)))
sum(y^2)
},method="BFGS")$par
return(list(f1=(((1/(del/sigma+(lambda*par[2]*del/sigma+del)/(par[6]*m*(par[3]^(alpha)))+(del+par[7]*m*(par[3]^(alpha)))/(lambda*par[1])+1)))*((lambda*par[2]*del/sigma+del)/(par[6]*m*(par[3]^(alpha)))))*((beta*0.5*v*(1-par[1])/(r+lambda)+((sister-(eta*(par[6]^t)/t))+par[6]*m*(par[3]^(alpha))*beta*(1-par[2])/(r+lambda))/r+(baby-sister-(eta*(par[7]^t-par[6]^t)/t)+par[7]*m*(par[3]^(alpha))*(beta*0.5*v*(1-par[1])/(r+lambda))-par[6]*m*(par[3]^(alpha))*(beta*0.5*v*(1-par[2])/(r+lambda)))/(r+del)))+(((1/(del/sigma+(lambda*par[2]*del/sigma+del)/(par[6]*m*(par[3]^(alpha)))+(del+par[7]*m*(par[3]^(alpha)))/(lambda*par[1])+1)))*del/sigma)*((beta*0.5*v*(1-par[2])/(r+lambda)+((sister-(eta*(par[6]^t)/t))+par[6]*m*(par[3]^(alpha))*beta*(1-par[2])/(r+lambda))/r))+((1/(del/sigma+(lambda*par[2]*del/sigma+del)/(par[6]*m*(par[3]^(alpha)))+(del+par[7]*m*(par[3]^(alpha)))/(lambda*par[1])+1)))*((((sister-(eta*(par[6]^t)/t))+par[6]*m*(par[3]^(alpha))*beta*(1-par[2])/(r+lambda))/r+(baby-sister-(eta*(par[7]^t-par[6]^t)/t)+par[7]*m*(par[3]^(alpha))*(beta*0.5*v*(1-par[1])/(r+lambda))-par[6]*m*(par[3]^(alpha))*(beta*0.5*v*(1-par[2])/(r+lambda)))/(r+del)))+(((lambda*par[2]*del/sigma+del)/(par[6]*m*(par[3]^(alpha))))*((1/(del/sigma+(lambda*par[2]*del/sigma+del)/(par[6]*m*(par[3]^(alpha)))+(del+par[7]*m*(par[3]^(alpha)))/(lambda*par[1])+1))))*((((sister-(eta*(par[6]^t)/t))+par[6]*m*(par[3]^(alpha))*beta*(1-par[2])/(r+lambda))/r))
))
}
sigma=seq(0.1,0.8,length.out=8)
del=seq(0.1,0.8,length.out=8)
baby=seq(0.1,0.8,length.out=8)
sister=seq(0.1,0.8,length.out=8)
z=array(0,c(1,8,8,8,8))
for(i in 1:8)
{for(j in 1:8)
{for(n in 1:8)
{for(l in 1:8)
{
gg=g(sigma[i],del[j],baby[n],sister[l])
z[1,i,j,n,l]=gg$f1
}
}
}
}