% The script file pwsemb04 calculates examples of once and
% twice differentiable processes with no dominant frequency.

dt=0.01; eta=1/5; n=300; etadt=eta*dt; gamma=exp(-etadt); 
FI=gamma*[1,0;etadt,1]; 
sigma1=1; sigma0=sigma1*sqrt(2);
P=[sigma0*sigma0,sigma1*sigma1;sigma1*sigma1,sigma1*sigma1];
p11=sqrt(P(1,1)); p12=0; p21=P(1,2)/p11; p22=sqrt(P(2,2)-p21*p21);
p=[p11,p12;p21,p22];
QQ=P-FI*P*FI';
q11=sqrt(QQ(1,1)); q21=QQ(1,2)/q11; q12=0; q22=sqrt(QQ(2,2)-q21*q21);
Q=[q11,q12;q21,q22]; 
ap=p*[randn;randn]; MA=[ap'];
   
for r=1:n-1
   ab=FI*ap+Q*[randn;randn];
   MA=[MA;ab'];
   ap=ab;
end      
      
   
T=0:dt:(n-1)*dt;T=T';
   
plot(T,MA(:,1),'k',T,MA(:,2),'r')
title('The once differentiable function and the auxillary one')
pause
close

dt=0.01; eta=1/5; n=300; etadt=eta*dt; gamma=exp(-etadt); 

FI=gamma*[1,0,0;etadt,1,0;etadt*etadt/2,etadt,1]; 
sigma2=1; 
P=sigma2*sigma2*[8/3,4/3,2/3;4/3,4/3,1;2/3,1,1];
p11=sqrt(P(1,1));p12=0;p13=0;p21=P(1,2)/p11;p31=P(1,3)/p11; 
p22=sqrt(P(2,2)-p21*p21);p23=0; p32=(P(2,3)-p31*p21)/p22;
p33=sqrt(P(3,3)-p31*p31-p32*p32);
p=[p11,p12,p13;p21,p22,p23;p31,p32,p33];
QQ=P-FI*P*FI';
q11=sqrt(QQ(1,1)); q21=QQ(1,2)/q11; q31=QQ(1,3)/q11; q12=0; q13=0;
q22=sqrt(QQ(2,2)-q21*q21);q23=0; q32=(QQ(2,3)-q31*q21)/q22;
q33=sqrt(QQ(3,3)-q31*q31-q32*q32);

Q=[q11,q12,q13;q21,q22,q23;q31,q32,q33]; 
ap=p*[randn;randn;randn]; MA=[ap'];
   
for r=1:n-1
   ab=FI*ap+Q*[randn;randn;randn];
   MA=[MA;ab'];
   ap=ab;
end      
      
   
T=0:dt:(n-1)*dt;T=T';
plot(T,MA(:,2),'k',T,MA(:,3),'r')
title('The twice differentiable function and the auxillary one')
pause
close

