Help needed how to let the code run on jacket and speed up!

[Old posts from the commercial version of ArrayFire] Discussion of ArrayFire using CUDA or OpenCL.

Moderator: pavanky

Help needed how to let the code run on jacket and speed up!

Postby henry.lou » Sun Aug 08, 2010 9:26 pm

Dear all:

I have a code and want let run on jacket. Where to modify and let the code speed up?
The code is below:

clear

tic
t=1
s=t+1 %賣掉時的日期,在Duffie and Lando (2000),s是當到期日
T=s+1 %當到期日



qq=100
qx=100
NN=100% 31230.219000 seconds seconds
Ny0=100
nox=100
dt=0.1
Q=10
d=0.05
aa3=[7.5+d:d:10]
z0=[log(80),log(90),log(99)]
bt=0.9
et=bt*qq
gy=qq-et

opport=[0.432414,0.452013,0.115573]%have changed% Rockafellar and Uryasev(1999):Gov bond低風險,投資比例低
numcorp=numel(z0)
z00=repmat(z0,qq,1)
%yy=[log(79);log(80);log(81);log(82);log(83);log(84);log(85);log(86);log(87);log(88);log(89);log(90);log(91);log(92);log(93);log(94);log(95);log(96);log(97);log(98);log(99);log(100)]


e4=log(78)
syms x

xx=linspace(e4,log(100),qx)'
xm=xx
se=0.05
var=se^2
m=0.01
k=se*t^0.5
noy=qq
%zzt=(z00+m*t+se*ran*t^0.5)
xi=repmat(x,1,numcorp)
%xm=log(linspace(1,300,nox))'
r=0.06
var=se^2
disc=exp(-r*(s-t))
%p=zeros(qq,numcorp)

numa=numel(aa3)


for ai=1:numa
a=aa3(1,ai);


e5=-0.5*a^2;


format long g

for yy=1:2
if yy==1
y=z0;

else

y=z0+m*t+0.5*var*t;%在t時的report y
end



for i=1:nox

syms x
nn=vpa((1-exp(-2*(z0-e4).*(xi-e4)./((se*t^0.5)^2))).*(normpdf(y-xi,e5,a)).*(normpdf(xi,m*t+z0,se*t^0.5)),5);
% when t is very small,1-exp(-2*(z0-e4)*(x-e4)/((se*t^0.5)^2))will be negative and not satisfy probablity axiom
ny=vpa((normpdf(y,m*t+z0+e5,(a^2+var*t)^0.5)),5);
b=vpa(nn./ny,5);

for bi=1:numcorp
pt(1,bi)=sum(subs(b(:,bi),x,xx))*((log(100)-e4)/qx)+0.00000000001;

end

x=xm(i,1);
% if 1-exp(-2*(z0-e4)*(x-e4)/((se*t^0.5)^2+0.1))<=0
% b=0
% end
gf(i,1:numcorp)=max(eval(vpa(b./pt,5)),0)+0.0000000000001;

end%i

totalgf=sum(gf);

for proi=1:numcorp
prox(:,proi)=gf(:,proi)./totalgf(1,proi);
end



for ii=1:Q
if yy==1
numx=prox.*Ny0;
fixx=fix(numx);
N=sum(fixx);
for coi=1:numcorp
newzz=[];
for c=1:nox

newzz=[newzz;repmat(xm(c,1),fixx(c,coi),1)];
end



for zi=1:N(1,coi)
zt=newzz(zi,1);

t1=0
nn=(s-t)/dt
if nn<=1
nn=1
end


for ni=1:nn

dw=randn*dt^0.5;
dz=m*dt+se*dw;
zt=(zt+dz);
zz=exp(zt);
s1=t1+dt;

if zz<=78
p1=0.4*exp(r*(t-s1));
%bre=bre+1
break
end
t1=t1+dt;
p1=1*exp(-r*(T-t));
end%ni

ppp(zi,coi)=vpa(p1,5);%不同y價格向量會有不同維度
end%zi
psta(ai,coi)=sum(ppp(1:N(1,coi),coi))/N(1,coi); %initial price

end%coi




else %yy=2



for j=1:qq

numx=prox.*NN;
fixx=fix(numx);
N=sum(fixx);
for coi=1:numcorp
newzz=[];
for c=1:nox

newzz=[newzz;repmat(xm(c,1),fixx(c,coi),1)];
end
for zi=1:N(1,coi)
zt=newzz(zi,1);
t1=t;
nn=(s-t)/dt
if nn<=1
nn=1
end

for ni=1:nn

dw=randn*dt^0.5;
dz=m*dt+se*dw;
zt=(zt+dz);
zz=exp(zt);
s1=t1+dt;

if zz<=78
p1=0.4*exp(r*(s-s1))
%bre=bre+1
break
end
t1=t1+dt;
p1=1*exp(-r*(T-s));
end%ni

pp(zi,1)=vpa(p1,5);%不同y價格向量會有不同維度
end%zi
avpp(ii,j,coi)=sum(pp(1:N(1,coi),1))/N(1,coi);
clear pp
end%coi

end %j

end %else if

end %ii


end%yy


% 由之前的結果可知,用模擬或模型公式,其估算出來的價格很接近,所以先用模擬的方式試試

format long g

for ii=1:Q
aavp=avpp(ii,:,:);
avp=[aavp(:,:,1)' aavp(:,:,2)' aavp(:,:,3)'];
ggk=sort(avp*opport');




%rankk=ggk(et,1)

af=opport*psta(ai,:)'-ggk(gy,1);


for k=1:qq
los(k,1)=vpa((psta(ai,:)-avp(k,1:numcorp))*opport',8);
lf(k,1)=max(0,eval(los(k,1)-af));
end

cvar=eval(af+1/(1-bt)*(1/qq)*sum(lf));
%load('CVaR0517a01','cvar01')

CVaR3(ii,ai)=cvar;
save('test0803c3a','CVaR3','a');
aff(ii,ai)=af;
clear af
end%ii

end %ai

cc3=mean(CVaR3)
save('test0803c3a','aff','avpp','y','z0','psta','CVaR3','cc3','aa3')

plot(aa3,cc3)
toc
henry.lou
 
Posts: 21
Joined: Thu Aug 13, 2009 6:18 am

Return to [archive-commercial] Programming & Development with ArrayFire

cron