b-s公式倒是明白，但是没用来计算过 credit risk，楼主的程序我给帮忙改了改，不知道是不是楼主的意思了

function[Va,SigmaVa]=CalV(E,SigmaE,D,r,T)

%a=zeros(8,8);

D=[

1096853709

1422513664

1080412570

1099330560

821720

676976

196840

37877619

];

a=D;

%b=zeros(8,8);

E=[

5482334

112522064

57679892

4537800

751520

1882103.375

285241.125

3833159

];

b=E;

%c=zeros(8,8);

SigmaE=[

0.401835787

0.543863954

0.609681884

0.741317746

0.917986928

0.931843334

1.004589468

1.02883818

];

c=SigmaE;

%计算Va'SigmaVa

key=0;

r=0.0225;

T=1;

PI=4*atan(1);

Va=b;%va以E为迭代的初值

newVa=Va;

SigmaVa=c; %SigmaVa以SigmaE为迭代的初值

for k=1:10  %迭代SigmaVa

for j=1:10 %迭代Va

oldVa=newVa;

% Va=(E+D*exp(-r,T)*Nd2)/Ndl;

% SigmaVa=E*SigmaE/(E+D*exp(-r,T)*Nd2);

%t1=log(oldVa/a);

d1=(log(oldVa./a)+(r+(SigmaVa.^2)/2)*T)./(SigmaVa.*(T^(1/2)));

d2=d1-(SigmaVa*(T^(1/2)));

Nd1=normcdf(d1,0,1);

Nd2=normcdf(d2,0,1);

fV=oldVa.*Nd1-a.*exp(-r*T).*Nd2-b;

temp111=exp(-d1.^2/2)./((2*PI*T)^0.5.*SigmaVa)

dfV=Nd1+temp111-a.*exp(-r*T).*exp(-d2.^2/2)./((2*PI*T)^0.5.*oldVa.*SigmaVa);

if dfV==0

dfV

pause;

end

dV=-fV./dfV;

newVa=oldVa+dV;   %前后两次的值认为非常接近，则迭代成功

newVa

if newVa==0

oldVa

newVa

pause;

end

if(dV/newVa)<1e-7 & (dV/newVa)>-1e-7

break;

end

end

Va=newVa;

d1=(log(Va./a)+(r+(SigmaVa.^2)/2)*T)./(SigmaVa.*(T^(1/2)));

d2=d1-(SigmaVa*(T^(1/2)));

Nd1=normcdf(d1,0,1);

Nd2=normcdf(d2,0,1);  %方程右边与左边的比值

diff=(Va.*Nd1-a.*exp(-r*T).*Nd2)./b;

if diff<1-1e-5 | diff>1+1e-5

assert('迭代出现错误'); %迭代出现错误

else

key=key+1;

end

newSigmaVa=b.*c./(b+a.*exp(-r*T).*Nd2);

if((SigmaVa-newSigmaVa)./newSigmaVa)<1e-5 & ((SigmaVa-newSigmaVa)./newSigmaVa)>(-1e-5)

break;

key=key+1;

else

SigmaVa=newSigmaVa;

end

end

if key==2

%else

% Va=一l：

% SigmaVa=-1：

end