Exercise - Generation of synthetic hydrological series for the Po River at Piacenza
This is a possible sequence of instructions in R. The list of instructions is uncommented.
It only serves as a guidance. The list assumes that the data are already uploaded in a
vector called: "po". Therefore, please note that the data need to be uploaded accordingly.


plot(po,type="l")
mean(po)
min(po)
max(po)
po1=log(po)
plot(po1,type="l")

po2=array(po1,dim=c(365,86))
mu=rep(0,365)
sigma=rep(0,365)

for(i in 1:365)
mu[i]=mean(po2[i,])
plot(mu,type="l")

for(i in 1:365)
sigma[i]=sd(po2[i,])
plot(sigma,type="l")

plot(sigma*10,type="l")
lines(mu,col="red")

po3=(po1-rep(mu,86))/rep(sigma,86)
plot(po3,type="l")
acf(po1,lag.max=365) # autocorrelation function useful to show the presence of seasonality
acf(po1,lag.max=365*87)

acf(po3,lag.max=365*87)
# plot(acf(po1),acf(po3), type="l",col="black","red")

phi1=acf(po3,lag.max=730)$acf[2]
phi1

phi1=acf(po3,lag.max=730)$acf[1]
phi1

phi1=acf(po3,lag.max=730)$acf[2]
epsilon=rep(0,365*86)
for(i in 2:(365*86))
epsilon[i]=po3[i]-phi1*po3[i-1]
plot(epsilon,type="l")
plot(density(epsilon),type="l")
acf(epsilon,lag.max=365)
acf(epsilon)

pr1=ar(po3,order.max=2) # change 2 moving forward 3,4,5....
names(pr1)
acf(pr1$resid[3:(365*86)]) # if changed above, increase also here 

mean(epsilon)
sd(epsilon)

epsilon1=rnorm(365*1000,0,sd(epsilon))
plot(epsilon1,type="l")

pos1=rep(0,(365*1000))
for(i in 2:(365*1000))
pos1[i]=phi1*pos1[i-1]+epsilon1[i]
plot(pos1[1:3650],type="l")
acf(pos1[1:3650],lag.max=365)

pos2=pos1*rep(sigma,1000)+rep(mu,1000)
acf(pos2[1:3650],lag.max=365)

pos3=exp(pos2)

mean(po)
mean(pos3)
sd(po)
sd(pos3)
max(po)
max(pos3)

plot(density(po),type="l")
lines(density(pos3),col="red",lwd=2) # peak lower due to distribution not gaussian

# pr1=acf(pos3)$acf # $acf for the autocorrelation coefficients
pr1=acf(pos3,lag.max=365)$acf
acf(po,lag.max=365)
lines(pr1,col="red",lwd=2) # for droughts is very fine