%%%%% Simuliamo l'equazione 
%%%%% dx(t) = -lambda x(t) dt + sigma dB(t)
%%%%% usando lo schema di Eulero esplicito:
%%%%% x[k+1]=(1-lambda*h)* x[k]+sqrt(h)* sigma* rnorm(1,0,1)

T<-100
h<- 0.1
lambda<- 1
sigma<- 1
x <- matrix(nrow=T/h)
x[1]<-0

for (k in 1:(T/h-1)) {
x[k+1]=(1-lambda*h)* x[k]+sqrt(h)* sigma* rnorm(1,0,1)

}
ts.plot(x)


%%%%% Supponiamo di avere la traccia "traccia1.txt",
%%%%% di cui ci viene dichiarato il tempo: 1 sec ogni 10 valori.
%%%%% Come scegliamo lambda e sigma?

x.emp<-read.table(file="traccia.txt")

ts.plot(x.emp)

hist(x.emp[,1],50)

qqnorm(x.emp[,1])

%%%%% è visibilmente gaussiana, il modello è ragionevole.
%%%%% Confrontiamo le acf:

acf(x.emp)

acf(x)

% Muoviamo per tentativi lambda per ottenere uno simile:

lambda<- 2
sigma<- 1
x <- matrix(nrow=T/h)
x[1]<-0

for (k in 1:(T/h-1)) {
x[k+1]=(1-lambda*h)* x[k]+sqrt(h)* sigma* rnorm(1,0,1)

}
acf(x)

%%%%% si vede che il metodo è impreciso. Aumentiamo la precisione: 

T<-1000
h<- 0.1
lambda<- 4
sigma<- 1
x <- matrix(nrow=T/h)
x[1]<-0

for (k in 1:(T/h-1)) {
x[k+1]=(1-lambda*h)* x[k]+sqrt(h)* sigma* rnorm(1,0,1)

}
acf(x,20)


%%%%% così è più precisa, ma come valutare la differenza?


T<-1000
h<- 0.1
lambda<- 5
sigma<- 1
x <- matrix(nrow=T/h)
x[1]<-0

for (k in 1:(T/h-1)) {
x[k+1]=(1-lambda*h)* x[k]+sqrt(h)* sigma* rnorm(1,0,1)

}
acf(x)[1:7]
acf(x.emp)[1:7]


%%%%% oppure usiamo la formula asintotica:

-(log(acf(x.emp)$acf))/((1:length(acf(x.emp)$acf))*h)


%%%%% da cui emerge un valore mediamente uguale a 5 se osserviamo
gli ultimi valori prima che si presentino gli zeri di acf.

lambda <- 5

%%%%% Stimiamo ora sigma, dalla formula:
%%%%% sigma^2/(2*lambda) = sd(x.emp[,1])^2

sigma <- sqrt(sd(x.emp[,1])^2*(2*lambda))
sigma

%%%%% Il valore vero era 3. Il fit non è precisissimo, ma c'è
l'ordine di grandezza.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Caso non gaussiano: densità esponenziale.
% Simuliamo l'equazione 
% dx(t) = (sigma-lambda*sigma/2*x(t))dt + sigma*sqrt(x(t))*dB(t)
% usando lo schema di Eulero esplicito:
% x[k+1]= x[k]+0.5*sigma^2*h 
% -0.5*lambda*sigma^2*x[k]*h
% +sigma*sqrt(h*x[k])*rnorm(1,0,1)
%%% Attenzione: per evitare i negativi, poniamo x[k+1]=0 se
%%% doveva essere negativo.

T<-10
h<- 0.01
lambda<- 1
sigma<- 1
xe <- matrix(nrow=T/h)
xe[1]<- 1
for (k in 1:(T/h-1)) {
A <-  xe[k]+0.5*sigma^2*h
B <-  -0.5*lambda*sigma^2*xe[k]*h
C <-  sigma*sqrt(h*xe[k])*rnorm(1,0,1)
xe[k+1]<-  (A+B+C)*(sign(A+B+C)+1)/2
}
ts.plot(xe)


%%% Controlliamo che la distribuzione al tempo 10 sia esponenziale:

T<-10
h<- 0.01
lambda<- 1
sigma<- 1

sample <- matrix(nrow=100)
for (i in 1:100) {

xe <- matrix(nrow=T/h)
xe[1]<- 1
for (k in 1:(T/h-1)) {
A <-  xe[k]+0.5*sigma^2*h
B <-  -0.5*lambda*sigma^2*xe[k]*h
C <-  sigma*sqrt(h*xe[k])*rnorm(1,0,1)
xe[k+1]<-  (A+B+C)*(sign(A+B+C)+1)/2
}

sample[i]<-xe[T/h]
}


hist(sample,20)


%%% Infine, potremmo scegliere il parametro sigma per "fittare"
%%% una certa acf empirica. Col sigma precedente vale:

acf(xe)

%%% Se il tasso di decadimento fosse lambda = 5, potremmo cercare 
%%% per tentativi:

T<-1000
h<- 0.01
lambda<- 1
sigma<- 3.2
xe <- matrix(nrow=T/h)
xe[1]<- 1
for (k in 1:(T/h-1)) {
A <-  xe[k]+0.5*sigma^2*h
B <-  -0.5*lambda*sigma^2*xe[k]*h
C <-  sigma*sqrt(h*xe[k])*rnorm(1,0,1)
xe[k+1]<-  (A+B+C)*(sign(A+B+C)+1)/2
}

-(log(acf(xe)$acf))/((1:length(acf(xe)$acf))*h)

%%% prendendo T più elevati la acf si stabilizza e la ricerca può
essere più precisa. 

T<-100
h<- 0.01
lambda<- 1
sigma<- 3.2
xe <- matrix(nrow=T/h)
xe[1]<- 1
for (k in 1:(T/h-1)) {
A <-  xe[k]+0.5*sigma^2*h
B <-  -0.5*lambda*sigma^2*xe[k]*h
C <-  sigma*sqrt(h*xe[k])*rnorm(1,0,1)
xe[k+1]<-  (A+B+C)*(sign(A+B+C)+1)/2
}

ts.plot(xe)