function [v,x,z,x1,w]=q4_8(n,b1)
% for simulations in question 4.8, b=2.3  and mu=-1.  
% run a bunch with    v = []; for i=1:50; v=[v;q4_8(1000,2.9)];   end
%  this function run and optimized using q4_8run
mu=-1; b=2.3;
t=trnd(3,5,n);
x=mu+t*(b/sqrt(3));
sumx=sum(x);
v=quantile(.95,sumx,ones(length(sumx)));
v=v(1);
%disp(['crude estimate of VaR ' num2str(v(1))])
%control variate estimate
u=tcdf(t,3);
z=norminv(u,mu,b);   % these are the normal control variates
v1=v;
for j=1:10
Gv=mean((sumx>v1)-(sum(z)>v1))+1-normcdf(v1,mu,b);   %estimates 1-F(v1)  using control variate.
v1=v1+(Gv-0.05)/normpdf(v1,mu,b); %  adjust the estimated VaR  using Fv
end
% importance sampling estimator. use as importance distribution a similar
% density but with parameters mu1,b1   
mu1=v/5;    x1=mu1+(b1/sqrt(3))*t;  % the importance variables
w=(b1^2 +(x1-mu1).^2)./(b^2 +(x1-mu).^2);
w=(b/b1)^15 *prod(w.^2);   % the vector of weights that attach to the sum(x1)
v2=v1;
for j=1:10
mn=weighted_mean((sum(x1)>v2),w);   %estimates 1-F(v2)  using control variate.
v2=v2+(mn-0.05)/normpdf(v2,mu1,b1); %  adjust the estimated VaR
end
%var(w.*(sum(x1)>v2))   % choose b1 (around 3) so that this is minimized
v=[v v1 v2];