%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% Exercice 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

mu=2;
x=0:.1:20;
lapdf=1/mu*exp(-x./mu);

% Graphe de la fonction cumulative
lacf=1-exp(-x./mu);
figure(1); clf;
subplot(211)
plot(x,lapdf)
subplot(212)
plot(x,lacf)

% Tirage de  N valeurs uniformement entre 0 et 1
N=10000;
val=rand(1,N);
va=-mu*log(1-val);

figure(2); clf;
subplot(211)
hist(va,100)

va2=exprnd(2,1,N);
subplot(212)
hist(va2,100)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% Exercice 2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all;

N=1000;
x=-20:.1:20;
sigma=1;
lapdf=1/sigma/sqrt(2*pi)*exp(-x.^2/2/sigma^2);


N=100;
Nreal=1000;
% Bruit blanc gaussien de N temps pour Nreal realisations
xxt=randn(Nreal,N);
mcovest=cov(xxt);
imagesc(mcovest)
% estimation fct de correlation (car stationnaire)
for ii=-N:N,
    cor(ii+N+1)=sum(diag(mcovest,ii))/(102-abs(ii));
end
figure;
plot(cor)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% Exercice 3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear ;
N=2^8;  	%100
Nreal=50;	%200
a=0.8;


%1% Programmation directe;
x=randn(Nreal,N);
y=zeros(Nreal,N); 
	% Note : etendre la taille d'une variable prend beaucoup
	% de temps a matlab ; il est mieux de lui donner la taille avant
y(:,1)=randn(Nreal,1)/sqrt(1-a^2);
	% initialisation en distribution Gaussienne qui sera 
	% la distribution stationnaire du processus
for ii=2:N,
y(:,ii)=a*y(:,ii-1)+x(:,ii);
end;

figure;
for jj=1:8,	%dessin de 8 realisations
subplot(8,1,jj);
plot(y(jj,:));
end;



%2-a% Programmation en resolvantl'equation par recurrence
x=randn(Nreal,N);
tic;
for j=1:Nreal
for ii=1:N
  y(j,ii)=sum(x(j,1:ii).*a.^(ii-1:-1:0));
end
end ; toc;
%2-b% Meilleure programmation en matlab de la meme chose:
vNreal=ones(Nreal,1);
y2=zeros(Nreal,N);
tic;
for ii=1:N,
y2(:,ii)=sum(x(:,1:ii).*(vNreal*a.^(ii-1:-1:0)),2);
end; toc;

% Estimation de la covariance
mcovest=cov(y);
imagesc(mcovest)
plot(y)

% Estimation fct correlation
for ii=-N:N,
    cor(ii+N+1)=sum(diag(mcovest,ii))/(N+1-abs(ii));
end
plot(cor)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% Exercice 4
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
close all
a=0.1
tau=0:1:100;

mcov=toeplitz(exp(-a.*tau.^2)); 
R=chol(mcov);

%% matrice de covariance
figure(1)
imagesc(mcov)

%%% decomposition de Cholesky
figure(2)
imagesc(R)

%% fn de correlation
plot(R(1,:))

%% generation de la va 
zz=randn(101,1);
xx=R'*zz;
plot(xx)

%% estimation de la matrice de covariance
clear xxt
for jj=1:100;
zz=randn(101,1);
xxt(:,jj)=R'*zz;
end;

xxt=xxt';
mcovest=cov(xxt);

close all
figure(1)
imagesc(mcovest)
figure(2)
imagesc(mcov)

%%% comparaison estime/prescrit
figure(3)
imagesc((mcov-mcovest).^2)

% estimation fct correlation
for ii=-101:101,
    cor(ii+102)=sum(diag(mcovest,ii))/(102-abs(ii));
end;

for ii=-101:101,
    cor2(ii+102)=sum(diag(mcov,ii))/(102-abs(ii));
end;
plot(-101:101,cor,'b',-101:101,cor2,'r')


% 
imagesc(min(mcovest,1))

%% Refaire la meme chose avec :
mcov=toeplitz(exp(log(a).*tau));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Exercice 5 : mouvement brownien 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
N=1000;

m_row=(1:N)'*ones(1,N);
m_col=ones(N,1)*(1:N);
mcov_mb=min(m_row,m_col);

Rmb=chol(mcov_mb);
zz=randn(N,1);
mb=Rmb'*zz;
plot(mb)

plot(1:N,cumsum(zz),1:N,mb,'.')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Exercice 6 : brownien fractionnaire
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
N=1000;
H=0.8;
H2=2*H;

m_row=((1:N)'*ones(1,N)).^H2;
m_col=(ones(N,1)*(1:N)).^H2;
toepfbm=toeplitz((0:N-1).^H2);

R=m_row+m_col-toepfbm;

Rchol=chol(R);
zz=randn(N,1);
mbf=Rchol'*zz;
figure;
plot(mbf) 

% moyenne et variance en chaque pt => non stationnaire

%% estimation pour le bgf
bgf=diff(mbf);

% Faire l'histogramme
% Evaluer moyenne et var en chaque pt (car processus stationnaire)
% Etude du processus stat
%   loi de proba
%   estimation fct de correlation
%


%%%%%%%%%%%%%%%%%%%%
%%% Exercice 7  %%%%
%%%%%%%%%%%%%%%%%%%%
%% A mettre en forme
%% effet de la convolution = filtrage

figure;
loglog(abs(fft(conv(mbf,1))))
hold
loglog(abs(fft(conv(mbf,gausswin(50)))))
loglog(abs(fft(conv(mbf,gausswin(10)))))
 
%%%%%%%%%%%%%%%%%%%%
%%% Exercice 8  %%%%
%%%%%%%%%%%%%%%%%%%%

% une realisation
N=2^16;
a=0.99; clear Snu;
ff=pi*(1:N)/N;
Snu(1:N/2)=(a./(a^2+4*pi^2*ff(1:N/2).^2));
Snu(N/2+1:N-1)=fliplr(Snu(1:N/2-1));
Snu(N)=0; Snu=Snu';

Zincr=fft(randn(N,1));
% ou bien (equivalent)
% Zincr=randn(N,1).*exp(i*2*pi*rand(N,1));
Xincr=Zincr.*sqrt(Snu);
xx=real(ifft(Xincr));

figure(1); plot(xx);
figure(2);
[SP,f]=pwelch(xx,100);
plot(f/2/pi,SP*3,'r',[0 ff(1:N/2)/pi],[Snu(N) Snu(1:N/2)'],'k');

N=1024;
M=500;
a=0.8;
ff=pi*(1:N)/N; 
clear Snu;
Snu(1:N/2)=(a./(a^2+4*pi^2*ff(1:N/2).^2));
Snu(N/2+1:N-1)=fliplr(Snu(1:N/2-1));
Snu(N)=0; Snu=Snu';
clear xxt
for jj=1:M;
    Zincr=fft(randn(N,1));
    Xincr=Zincr.*sqrt(Snu);
    xxt(:,jj)=real(ifft(Xincr));
end;

xxt=xxt';
mcovest=cov(xxt);
tau=0:N-1;
mcov=toeplitz((a).^tau)*mean(diag(mcovest));

figure(1); colorbar
imagesc(mcovest)
figure(2); colorbar
imagesc(mcov)
figure(3); colorbar
imagesc((mcov-mcovest).^2)

cor=[]; cor2=[];
for ii=-N/2:N/2,
    cor(ii+N/2+1)=sum(diag(mcovest,ii))/(N+1-abs(ii));
end;
for ii=-N/2:N/2,
    cor2(ii+N/2+1)=sum(diag(mcov,ii))/(N+1-abs(ii));
end;
figure(4); clf; 
plot(-N/2:N/2,cor,'b');
plot(-N/2:N/2,cor,'b',-N/2:N/2,cor2,'r');


