lines(GAcorr2@summary[,1],lty=9)
lines(GAcorr1@summary[,1],lty=10)
setwd("Desktop/Carla")
setwd("/Users/olmi/Desktop/Carla")
library(readxl)
read_xlsx("etichette.xlsx")
dati <- read_xlsx("etichette.xlsx")
names(dati)
head(dati)
GAcorr10@bestSol
GAcorr10@solution
ARMA.params <- function(crom){
mu <- crom[1]
p <- round(crom[2])
q <- round(crom[3])
fi <- c()
teta <- c()
if (p>0)
for (i in 1:p) fi = cbind(fi,crom[i+3])
if (q>0)
for (j in 1:q) teta = cbind(teta,crom[j+p+3])
print(mu)
print(c(p,fi))
print(c(1,teta))
}
ARMA_params(GAcorr10@solution)
ARMA.params(GAcorr10@solution)
ARMA.params <- function(crom){
mu <- crom[1]
p <- round(crom[2])
q <- round(crom[3])
fi <- c()
teta <- c()
if (p>0)
for (i in 1:p) fi = cbind(fi,crom[i+3])
if (q>0)
for (j in 1:q) teta = cbind(teta,crom[j+p+3])
print(mu)
print(c(p,fi))
print(c(q,teta))
}
ARMA.params(GAcorr10@solution)
ARMA.params <- function(crom){
mu <- crom[1]
p <- round(crom[2])
q <- round(crom[3])
fi <- c()
teta <- c()
if (p>0)
for (i in 1:p) fi = cbind(fi,crom[i+3])
if (q>0)
for (j in 1:q) teta = cbind(teta,crom[j+p+3])
print(c("mean",mu)
print(c(p,fi))
print(c(q,teta))
}
ARMA.params <- function(crom){
mu <- crom[1]
p <- round(crom[2])
q <- round(crom[3])
fi <- c()
teta <- c()
if (p>0)
for (i in 1:p) fi = cbind(fi,crom[i+3])
if (q>0)
for (j in 1:q) teta = cbind(teta,crom[j+p+3])
print(c("mean",mu))
print(c(p,fi))
print(c(q,teta))
}
ARMA.params(GAcorr10@solution)
ARMA.params <- function(crom){
mu <- crom[1]
p <- round(crom[2])
q <- round(crom[3])
fi <- c()
teta <- c()
if (p>0)
for (i in 1:p) fi = cbind(fi,crom[i+3])
if (q>0)
for (j in 1:q) teta = cbind(teta,crom[j+p+3])
print(c("mean",mu))
print(c("AR",p,fi))
print(c("MA",q,teta))
}
ARMA.params(GAcorr10@solution)
ARMA.params <- function(crom){
mu <- crom[1]
p <- round(crom[2])
q <- round(crom[3])
fi <- c()
teta <- c()
if (p>0)
for (i in 1:p) fi = cbind(fi,crom[i+3])
if (q>0)
for (j in 1:q) teta = cbind(teta,crom[j+p+3])
print(cat("mean",mu))
print(c("AR",p,fi))
print(c("MA",q,teta))
}
ARMA.params(GAcorr10@solution)
cat(1,2)
cat("mean",mu)
ARMA.params <- function(crom){
mu <- crom[1]
p <- round(crom[2])
q <- round(crom[3])
fi <- c()
teta <- c()
if (p>0)
for (i in 1:p) fi = cbind(fi,crom[i+3])
if (q>0)
for (j in 1:q) teta = cbind(teta,crom[j+p+3])
cat("mean",mu)
print(c("AR",p,fi))
print(c("MA",q,teta))
}
ARMA.params(GAcorr10@solution)
ARMA.params <- function(crom){
mu <- crom[1]
p <- round(crom[2])
q <- round(crom[3])
fi <- c()
teta <- c()
if (p>0)
for (i in 1:p) fi = cbind(fi,crom[i+3])
if (q>0)
for (j in 1:q) teta = cbind(teta,crom[j+p+3])
cat("mean",mu)
cat("AR",p,fi)
cat("MA",q,teta)
}
ARMA.params(GAcorr10@solution)
ARMA.params <- function(crom){
mu <- crom[1]
p <- round(crom[2])
q <- round(crom[3])
fi <- c()
teta <- c()
if (p>0)
for (i in 1:p) fi = cbind(fi,crom[i+3])
if (q>0)
for (j in 1:q) teta = cbind(teta,crom[j+p+3])
cat("mean",mu,"\n")
cat("AR",p,fi)
cat("MA",q,teta)
}
ARMA.params(GAcorr10@solution)
ARMA.params <- function(crom){
mu <- crom[1]
p <- round(crom[2])
q <- round(crom[3])
fi <- c()
teta <- c()
if (p>0)
for (i in 1:p) fi = cbind(fi,crom[i+3])
if (q>0)
for (j in 1:q) teta = cbind(teta,crom[j+p+3])
cat("mean",mu,"\n")
cat("AR",p,fi,"\n")
cat("MA",q,teta,"\n")
}
ARMA.params(GAcorr10@solution)
arima(X,order=c(1,0,1))
plot(GAcorr10)
help("plot.ga")
plot(GAcorr10,legend=FALSE)
setwd("/Users/olmi/Desktop/RPA-current/R code/ARMA")
setwd("/Users/olmi/Desktop/RPA-current/R code/ARMA")
save.image("~/Desktop/RPA-current/R code/ARMA/Calcolo_ARMA_Loire.RData")
plot(GAcorr10)
## Code_6_2.R
# Moving-average process
# MA(1) positive coefficient
# for reproducibiity
set.seed(13)
# generate white (Gaussian) noise
w <- rnorm(100)
x <- rep(0,100)
# generate MA(1)
x[1] <- w[1]
for (t in 2:100) x[t] <- w[t] + 0.5*w[t-1]
# plot time series
plot(as.ts(x),ylab="MA(1) process")
# plot autocorrelation function
acf(x,main="")
# plot spectral density
x.spec <- spectrum(x,log="no",span=30,plot=FALSE)
plot(x.spec$freq,x.spec$spec,t="l",ylim = c(0,3),
xlab="frequency (cycles/sample interval)",ylab="spectral density")
## Code_6_3.R
#Autoregressive model (1)
# Set seed for random number generation, for reproducibility of results
set.seed(45)
# Common parameters
N <- 101
t <- 0:(N-1)
# white (gaussian) noise with variance = 1
epsilon <- rnorm(N)
# AR(1) stationary
phi <- 0.3
Xs <- rep(0,N)
for (i in 2:N) {
Xs[i] <- phi*Xs[i-1] + epsilon[i]
}
# AR(1) not stationary
phi <- 1.01
X <- rep(0,N)
for (i in 2:N) {
X[i] <- phi*X[i-1] + epsilon[i]
}
# Comparison
plot(t,X,t="l")
lines(t,Xs,col="red")
## Code_6_4.R
#Autoregressive model (2) parameters
# AR(2)
phi1 <- 0.5
phi2 <- 0.2
X <- rep(0,N)
for (i in 3:N) {
X[i] <- phi1*X[i-1] + phi2*X[i-2] + epsilon[i]
}
plot(t,X,t="l")
## Code_6_5.R
#Autocorrelation structure of AR(2) and MA(2)
# stationary AR(2)
set.seed(4321)
x.AR2 <- arima.sim(model = list(order=c(2,0,0), ar = c(0.5,-0.7)),
n = 200)
plot(x.AR2,main="")
acf(x.AR2,main="")
pacf(x.AR2,main="")
# invertible MA(2)
set.seed(4321)
x.MA2 <- arima.sim(model = list(order=c(0,0,2), ma = c(0.2,-0.5)),
n = 200)
plot(x.MA2,main="")
acf(x.MA2,main="")
pacf(x.MA2,main="")
## Code_6_5.R
#Autocorrelation structure of AR(2) and MA(2)
# stationary AR(2)
set.seed(4321)
x.AR2 <- arima.sim(model = list(order=c(2,0,0), ar = c(0.5,-0.7)),
n = 200)
plot(x.AR2,main="")
acf(x.AR2,main="")
pacf(x.AR2,main="")
# invertible MA(2)
set.seed(4321)
x.MA2 <- arima.sim(model = list(order=c(0,0,2), ma = c(0.2,-0.5)),
n = 200)
plot(x.MA2,main="")
acf(x.MA2,main="")
pacf(x.MA2,main="")
## Code_6_8.R
# Difference operator
# for repeatability
set.seed(1234)
# generate data from an ARIMA(2,1,1) model
# AR parameters: 0.5, -0.7
# MA parameter: 0.6
x.ARIMA <- arima.sim(model = list(order=c(2,1,1), ar = c(0.5,-0.7), ma = 0.6), n = 200)
plot(x.ARIMA,main="",ylab="time series")
acf(x.ARIMA,main="")
x.ARMA = diff(x.ARIMA)
plot(x.ARMA,main="",ylab="time series")
acf(x.ARMA,main="")
aic.fit <- c()
for (p in 1:3){
for (q in 1:3) {
fit <- arima(x.ARMA,order = c(p,0,q))
aic.fit <- c(aic.fit,fit$aic)
}
}
AIC.fit <- matrix(data=aic.fit,nrow=3,ncol=3)
AIC.fit
arima(x.ARMA,order=c(2,0,1))
## Code_6_7.R
# Seasonal ARIMA
setwd("C:\RPA\code\ARMA\data")
# Maximum temperature in Bologna (Period: 1961-1991)
Tmax.data <- read.table("Bologna_Tmax_1961_1991.txt",h=FALSE)
Tmax <- Tmax.data$V1
# Compute lag-12 differences
Tmax.diff12 <- diff(Tmax,12)
# Plots
ts.plot(Tmax,xlab="Months from January 1961")
ts.plot(Tmax[1:60],xlab="Years 1961-1991",ylab="Tmax")
# ACF
acf(Tmax,lag.max = 36,main="")
acf(Tmax.diff12,lag.max = 36,main="")
pacf(Tmax.diff12,lag.max = 36,main="")
setwd("C:/RPA_FINAL/code/ARMA")
## Code_6_9.R
# Seasonal ARIMA
setwd("C:\RPA\code\ARMA\data")
# Maximum temperature in Bologna (Period: 1961-1991)
Tmax.data <- read.table("Bologna_Tmax_1961_1991.txt",h=FALSE)
Tmax <- Tmax.data$V1
# Compute lag-12 differences
Tmax.diff12 <- diff(Tmax,12)
# Plots
ts.plot(Tmax,xlab="Months from January 1961")
ts.plot(Tmax[1:60],xlab="Years 1961-1991",ylab="Tmax")
# ACF
acf(Tmax,lag.max = 36,main="")
acf(Tmax.diff12,lag.max = 36,main="")
pacf(Tmax.diff12,lag.max = 36,main="")
## Code_6_9.R
# Seasonal ARIMA
setwd("C:\RPA\code\ARMA\data")
## Code_6_9.R
# Seasonal ARIMA
setwd("C:RPA\code\ARMA\data")
setwd("C:/RPA/code/ARMA")
setwd("C:/RPA/code/ARMA")
# Maximum temperature in Bologna (Period: 1961-1991)
Tmax.data <- read.table("Bologna_Tmax_1961_1991.txt",h=FALSE)
Tmax <- Tmax.data$V1
# Compute lag-12 differences
Tmax.diff12 <- diff(Tmax,12)
# Plots
ts.plot(Tmax,xlab="Months from January 1961")
ts.plot(Tmax[1:60],xlab="Years 1961-1991",ylab="Tmax")
# ACF
acf(Tmax,lag.max = 36,main="")
acf(Tmax.diff12,lag.max = 36,main="")
pacf(Tmax.diff12,lag.max = 36,main="")
# Maximum temperature in Bologna (Period: 1961-1991)
Tmax.data <- read.table("Bologna_Tmax_1961_1991.txt",h=FALSE)
# Maximum temperature in Bologna (Period: 1961-1991)
Tmax.data <- read.table("data/Bologna_Tmax_1961_1991.txt",h=FALSE)
## Code_6_9.R
# Seasonal ARIMA
setwd("C:/RPA/code/ARMA")
# Maximum temperature in Bologna (Period: 1961-1991)
Tmax.data <- read.table("data/Bologna_Tmax_1961_1991.txt",h=FALSE)
Tmax <- Tmax.data$V1
# Compute lag-12 differences
Tmax.diff12 <- diff(Tmax,12)
# Plots
ts.plot(Tmax,xlab="Months from January 1961")
ts.plot(Tmax[1:60],xlab="Years 1961-1991",ylab="Tmax")
# ACF
acf(Tmax,lag.max = 36,main="")
acf(Tmax.diff12,lag.max = 36,main="")
pacf(Tmax.diff12,lag.max = 36,main="")
## Code_6_10.R
# best fit SARIMA orders
p = 1
q = 0
P = 1
Q = 1
# SARIMA fit
(arima(Tmax,order=c(p,0,q),list(order=c(P,1,Q),period=12)) -> fit)
acf(residuals(fit),48,main="")
hist(residuals(fit),main="")
## Code_6_11.R
# convert Tmax into a time series
# X.true is the full time series
X.true <- ts(Tmax)
# X is the time series of the period 1961-1989
X <- ts(Tmax[1:(29*12)],start=1,end=(1989-1961+1)*12)
# seasonal ARIMA model
fit <- arima(X,order=c(1,0,0),list(order=c(1,1,1),period=12))
# make prediction for next 2 years
X.pred <- predict(fit, 24)
# Plot predicted values
plot(X[313:348],t="l",xlim=c(1,60),xlab="January 1987 to December 1991",ylab="Tmax")
lines(37:60,X.pred$pred,col="red")
lines(36:60,X.true[(313+35):372],lty=2)
setwd("C:/RPA/code/ARMA")
## Code_6_12.R
#Loire river analysis
setwd("C:/RPA/code/ARMA")
read.table("data/Loire_rainfall1.txt",h=T) -> rain1
read.table("data/Loire_rainfall2.txt",h=T) -> rain2
# annual data
plot(rain1$Y,rain1$TOT,t="l",xlab="Year",ylab="Precipitation (mm)")
lines(rain2$Y,rain2$TOT,col="red")
# average between locations
rain.AVG <- (rain1$TOT+rain2$TOT)/2
hist(rain.AVG,freq=FALSE,main="",xlab="Total rain (average between locations)")
x <- rain.AVG
curve(dnorm(x, mean=mean(rain.AVG), sd=sd(rain.AVG)), add=TRUE, col="red")
qqnorm(x,main="")
qqline(x,col="red")
setwd("C:/RPA/code/ARMA")
## Code_6_12.R
setwd("C:/RPA/code/ARMA")
# read data
dati <- read.table("data/loire_runoff.txt",h=TRUE)
# convert discharge data to time series
X <- ts(dati$DISCHRG)
plot(X,ylab=expression(paste("Monthly runoff (m"^"3","/s)")),
xlab="Month # from January 1863",las=3)
# seasonal ARIMA analisys on monthly data
arima(X,order=c(1,0,2),seasonal=list(order=c(2,1,1),period=12)) -> fit
X.fitted <- X-fit$residuals
acf(X,na.action=na.pass,44,main="")
pacf(X,na.action=na.pass,44,main="")
plot(X.fitted,fit$residuals,xlab="Fitted",ylab="Residuals")
# seasonal ARIMA analisys on log of monthly data
plot(log(X),ylab=expression(paste("log(Monthly runoff) (m"^"3","/s)")),
xlab="Month # from January 1863",las=3)
arima(log(X),order=c(1,0,2),seasonal=list(order=c(2,1,1),period=12)) -> fit.log
lines(fitted(fit.log),col="red")
plot(fitted(fit.log),fit.log$residuals,xlab="Fitted",ylab="Residuals")
setwd("C:/RPA/code/ARMA")
# read data
dati <- read.table("data/loire_runoff.txt",h=TRUE)
# convert discharge data to time series
X <- ts(dati$DISCHRG)
plot(X,ylab=expression(paste("Monthly runoff (m"^"3","/s)")),
xlab="Month # from January 1863",las=3)
# seasonal ARIMA analisys on monthly data
arima(X,order=c(1,0,2),seasonal=list(order=c(2,1,1),period=12)) -> fit
X.fitted <- X-fit$residuals
acf(X,na.action=na.pass,44,main="")
pacf(X,na.action=na.pass,44,main="")
plot(X.fitted,fit$residuals,xlab="Fitted",ylab="Residuals")
# seasonal ARIMA analisys on log of monthly data
plot(log(X),ylab=expression(paste("log(Monthly runoff) (m"^"3","/s)")),
xlab="Month # from January 1863",las=3)
arima(log(X),order=c(1,0,2),seasonal=list(order=c(2,1,1),period=12)) -> fit.log
lines(fitted(fit.log),col="red")
plot(fitted(fit.log),fit.log$residuals,xlab="Fitted",ylab="Residuals")
## Code_6_12.R
setwd("C:/RPA/code/ARMA")
# read data
dati <- read.table("data/loire_runoff.txt",h=TRUE)
# convert discharge data to time series
X <- ts(dati$DISCHRG)
plot(X,ylab=expression(paste("Monthly runoff (m"^"3","/s)")),
xlab="Month # from January 1863",las=3)
# seasonal ARIMA analisys on monthly data
arima(X,order=c(1,0,2),seasonal=list(order=c(2,1,1),period=12)) -> fit
X.fitted <- X-fit$residuals
acf(X,na.action=na.pass,44,main="")
pacf(X,na.action=na.pass,44,main="")
plot(X.fitted,fit$residuals,xlab="Fitted",ylab="Residuals")
# seasonal ARIMA analisys on log of monthly data
plot(log(X),ylab=expression(paste("log(Monthly runoff) (m"^"3","/s)")),
xlab="Month # from January 1863",las=3)
arima(log(X),order=c(1,0,2),seasonal=list(order=c(2,1,1),period=12)) -> fit.log
lines(fitted(fit.log),col="red")
plot(fitted(fit.log),fit.log$residuals,xlab="Fitted",ylab="Residuals")
fit.log
lines(fitted(fit.log),col="red")
fitted(fit.log)
## Code_6_12.R
setwd("C:/RPA/code/ARMA")
# read data
dati <- read.table("data/loire_runoff.txt",h=TRUE)
# convert discharge data to time series
X <- ts(dati$DISCHRG)
plot(X,ylab=expression(paste("Monthly runoff (m"^"3","/s)")),
xlab="Month # from January 1863",las=3)
# seasonal ARIMA analisys on monthly data
arima(X,order=c(1,0,2),seasonal=list(order=c(2,1,1),period=12)) -> fit
X.fitted <- X-fit$residuals
acf(X,na.action=na.pass,44,main="")
pacf(X,na.action=na.pass,44,main="")
plot(X.fitted,fit$residuals,xlab="Fitted",ylab="Residuals")
# seasonal ARIMA analisys on log of monthly data
plot(log(X),ylab=expression(paste("log(Monthly runoff) (m"^"3","/s)")),
xlab="Month # from January 1863",las=3)
arima(log(X),order=c(1,0,2),seasonal=list(order=c(2,1,1),period=12)) -> fit.log
lines(fitted(fit.log),col="red")
plot(fitted(fit.log),fit.log$residuals,xlab="Fitted",ylab="Residuals")
## Code_6_12.R
setwd("C:/RPA/code/ARMA")
# read data
dati <- read.table("data/loire_runoff.txt",h=TRUE)
# convert discharge data to time series
X <- ts(dati$DISCHRG)
plot(X,ylab=expression(paste("Monthly runoff (m"^"3","/s)")),
xlab="Month # from January 1863",las=3)
# seasonal ARIMA analisys on monthly data
arima(X,order=c(1,0,2),seasonal=list(order=c(2,1,1),period=12)) -> fit
X.fitted <- X-fit$residuals
acf(X,na.action=na.pass,44,main="")
pacf(X,na.action=na.pass,44,main="")
plot(X.fitted,fit$residuals,xlab="Fitted",ylab="Residuals")
# seasonal ARIMA analisys on log of monthly data
plot(log(X),ylab=expression(paste("log(Monthly runoff) (m"^"3","/s)")),
xlab="Month # from January 1863",las=3)
arima(log(X),order=c(1,0,2),seasonal=list(order=c(2,1,1),period=12)) -> fit.log
lines(fitted(fit.log),col="red")
plot(fitted(fit.log),fit.log$residuals,xlab="Fitted",ylab="Residuals")
## Code_6_12.R
setwd("C:/RPA/code/ARMA")
# read data
dati <- read.table("data/loire_runoff.txt",h=TRUE)
# convert discharge data to time series
X <- ts(dati$DISCHRG)
plot(X,ylab=expression(paste("Monthly runoff (m"^"3","/s)")),
xlab="Month # from January 1863",las=3)
# seasonal ARIMA analisys on monthly data
arima(X,order=c(1,0,2),seasonal=list(order=c(2,1,1),period=12)) -> fit
X.fitted <- X-fit$residuals
acf(X,na.action=na.pass,44,main="")
pacf(X,na.action=na.pass,44,main="")
plot(X.fitted,fit$residuals,xlab="Fitted",ylab="Residuals")
# seasonal ARIMA analisys on log of monthly data
plot(log(X),ylab=expression(paste("log(Monthly runoff) (m"^"3","/s)")),
xlab="Month # from January 1863",las=3)
arima(log(X),order=c(1,0,2),seasonal=list(order=c(2,1,1),period=12)) -> fit.log
lines(fitted(fit.log),col="red")
plot(fitted(fit.log),fit.log$residuals,xlab="Fitted",ylab="Residuals")
