function [t0,p0]=DAmix(x)
% Data augmentation algorithm (see Liu p 134 ) applied to a normal variance mixture with
% 2 components
% test using load missing2; [t,p]=DAmix(z);
%e.g. [t,p1]=DAmix(normmixrnd(50,[0 0],[1 3],[2/3 1/3]));
t=[.000072 .000225];  % preliminary estimates of the variance of the two components
p1=.81;   t0=[]; p0=[];% preliminary estimates of the two probabilities
f1=normpdf(x,0,sqrt(t(1))); f2=normpdf(x,0,sqrt(t(2)));
w=(f1+(1/p1-1)*f2);
w=f1./max(w,.000001);
delta=(abs(x)<median(abs(x)));  n=length(x); ss=sum(x.^2);    
for i=1:10000
    n1=sum(delta); ss1=sum(delta.*x.^2); n2=sum(1-delta); ss2=sum((1-delta).*x.^2);
    if((ss1/n1)<(ss2/n2))   % do not run unless parameters are ok
     beta=2./max(.00001,[ss1 ss2]);
    %n1*beta(1)
y=gamrnd([n1 n2]/2,beta); t=1./y;
%p1=betarnd(n1+.5,n2+.5);  % percent this out if we want to fix p1.
f1=normpdf(x,0,sqrt(t(1))); f2=normpdf(x,0,sqrt(t(2)));
w=(f1+(1/p1-1)*f2);
w=f1./max(w,.000001);
delta=binornd(1,w);
%sum(delta)
t0=[t0;t]; p0=[p0 p1];   %    ss0=[ss0;[ss1/n1 ss2/n2]];
end
end
%plot annual volatility
 [N,X]=hist(sqrt(252*t0),30);
 plot(X,N)
 %figure
 %hist(p0,30)