function [p,beta,psi,ps]=stableexample(x,theta,a,m)
% x the data is  a vector of length 39, assumed to be stable(0,1,a) and we
% wish to estimate a-input an initial estimator
% e.g. [p,beta,psi,ps]=stableexample(stabrnd(1.5,39),1,1.5,10000)   m=number of
% simulations. 
ps=[]; de=[];
psi=[sum(trim(x-theta,.5)) sum(trim(x-theta,1)) sum(trim(x-theta,1.5)) sum(trim(x-theta,2))];
for i=1:m  % use m simulations to center and obtain covariance matrix etc.
z = theta+stabrnd(a,39);
est=[sum(trim(z-theta,.5)) sum(trim(z-theta,1)) sum(trim(z-theta,1.5)) sum(trim(z-theta,2))];
ps=[ps;mean(est')];
delt=[sum((abs(z-theta)<.5)) sum((abs(z-theta)<1)) sum((abs(z-theta)<1.5)) sum((abs(z-theta)<2))];
de=[de;delt];
end
C=cov(ps);
de=mean(de);
beta=de*inv(C);
ps=mean(ps);
p=beta*(psi-ps)';