function [par,op,param,SS]=calibratenig2(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]  ;
% same as calibratnig only we use linear regression
% [param,SS]=calibratenig2([4.4559    1.9958    0.8211])
% these parameters correpond to op=[168.6931  118.7571   83.0683   51.2564 38.8448   28.7075   20.7076 14.5844   10.0209    6.7354]

p=length(par); 
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
ns=p+1; np=p+1;  % start out with as many simulations ns as parameters np
param=repmat(par,ns, 1);  nsim=30000; seed=1;
param=param.*(.9+.2*rand(ns,p));   Y=zeros(ns,no);   % Y is ns by no  and i'th row is i'th parameter value
for i=1:ns   % 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
% now find a candidate z for the optimum
 SS=(Y-repmat(ca,np,1)).^2;
[s,k]=min(sum(SS,2));
 tn=param(k,:);   % this is the best parameter used so far-expand about this one
  stepsize=.1*sqrt(sum(tn.^2));  % sets the stepsize equal to 10% of the length of this vector
  %initialize
    a=zeros(ns,no); b=zeros(p,ns,no); grad=zeros(p,no); 
    SS=sum(SS,2);   % this is the vector of the sum of squared errors for each simulation.
while (size(Y,1)<50)
    for i=1:ns  % from here on repeat until enough simulations are done
    for j=1:no    % for each parameter param(i,:) value and each option j we estimate a,b
    [a(i,j),b(:,i,j),tk]=linreg(Y(:,j),param,param(i,:));  % fits a linear function for the jth option
end
end
% estimate the function at the new parameter value tn.  Fix the value of tk
Z=repmat(tn,ns,1); w=sum((param-Z).^2,2);  w=1./(1+tk*w); w=w/sum(w); % w is the weight vector
    % estimate the values of Pj  at the parameter value tn
    Phat=[];
    for j=1:no
        C=(Z-param)'.*b(:,:,j)*w;
        Phat=[Phat w'*a(:,j)+sum(C)];
    end
    derPhat=zeros(p,no);
    for j=1:no
    for k=1:ns    
        derPhat(:,j)=derPhat(:,j)+w(k)*b(:,k,j);
    end
end
    gradhat=2*derPhat*(Phat-ca)';
    gradhat=gradhat/sqrt(sum(gradhat.^2));
    tn=tn-stepsize*gradhat';   % the new parameter value 
    % add the next batch of simulations
    param=[param;tn];      ns=ns+1;    
par=max(0,tn);    % make sure parameters are positive
    alpha=abs(par(2)+1)+exp(par(1));
    op=nigcall(S0,K,.01,T,alpha,par(2),par(3),nsim,seed);  % the new option price corresponding to tn
    Y=[Y;op];     
    nsim=nsim+2000;   % increase the number of simulations at each step 
    SS=[SS ;sum((op-ca).^2)];
    if SS(ns)<SS(ns-1)
        stepsize=stepsize*(1-stepsize);
    end
    seed=seed+1000;
end
[s,k]=min(SS); 
par=param(k,:);
    alpha=abs(par(2)+1)+exp(par(1));
    op=nigcall(S0,K,.01,T,alpha,par(2),par(3),nsim,seed);