function [tn,SS]=calibratenig(par)
% uses calibration scheme in chapter 8
% based on the sum of squared distances to the observed options prices for NIG
% model. since we need alpha>|beta+1|  we write par=[log(alpha-abs(1+beta), beta,delta]  ;
% run this as....
% par=initial guess at parameter values. 
p=length(par);  % the number of parameters in a quadratic model is (p^2+3*p+2)/2 


ca=[173.1000  124.8000   88.5000   53.7000   39.3000   27.3000 17.8500 10.9500  6.2000   3.2500]; % observed prices
K=[950  1005 1050  1100  1125  1150   1175 1200  1225 1250]; % strike prices
no=length(K);  % no=number of options.
T=.3699; S0=1116.84;   % maturity and initial value
np=(p^2+3*p+2)/2;   % start out with as many simulations as parameters np
param=repmat(par,np, 1);  nsim=20000; seed=1;
param=param.*(.9+.2*rand(np,p));   Y=zeros(np,no);
for i=1:np   % index over the different options 
    par=param(i,:);
    alpha=abs(par(2)+1)+exp(par(1));
    op=nigcall(S0,K,.01,T,alpha,par(2),par(3),nsim,seed);
    Y(i,:)=op;  % 'th row of Y is vector of prices of the options for the ith value of the parameter
end


while (size(Y,1)<40)
    % now find a candidate z for the optimum
SS=(Y-repmat(ca,np,1)).^2;
[s,k]=min(sum(SS,2));
 z=param(k,:);   % this is the best parameter used so far-expand about this one
g=[];  CT=zeros(p,p); CBT=zeros(p,1); tn=z;
for l=1:2  % 2= number of times we update estimate of theta
for j=1:no
    [a,b,C,m]=quadreg(Y(:,j),param,z);  % fits a quadratic for the jth option
    eps=a+(tn-z)*b+0.5*(tn-z)*C*(tn-z)'-ca(j);
    CT=CT+eps*C; CBT=CBT+eps*b;
end
tn=max(0,real(z-CBT'*inv(CT))) ; % this is the parameter value for the next round of simulations.
end
% add the next batch of simulations
param=[param;tn];      np=np+1;    
par=max(0,param(np,:));
    alpha=abs(par(2)+1)+exp(par(1));
    op=nigcall(S0,K,.01,T,alpha,par(2),par(3),nsim,seed);
    Y=[Y;op]; 
    nsim=nsim+2000;
    sum((op-ca).^2)
    seed=seed+1000;
end
SS=sum(SS,2);