function [y,s]=hwvol(y0,rho,xi,k,ybar)
%simulates the squared volatiltiy in the Hull White model as specicied in Reesor&Mcleish
% y0= 0.1414; rho = -.1106; xi= 1.3671; k = 3.081,ybar=.1;
[w,s]=simbm(1,[0 1]); delta=1/10000;
phi=exp(-(k+.5*xi.^2).*s+xi*w);
y=y0*phi+k*ybar*cumsum(1./phi)*delta;
plot(s,y)
x=.01:.01:5
gx=x^(k/xi^2).*exp(