POSIXct (calendar time) e POSIXlt (local time)

# ISOdate(year, month, day, hour = 12, min = 0, sec = 0, tz = "GMT")

d <- ISOdate(2006,6,9)
d # [1] "2006-06-09 12:00:00 GMT"
class(d) # "POSIXt"  "POSIXct"

ISOdatetime(2006,6,9) # errore: l'ora non e' impostata come default a 12:00:00


as.POSIXlt(d) # in questa caso coincide
file.info(dir()) 
file.info(dir())[,"mtime"] # il concetto di time zone e' importante!

format

format(d,"%a") # giorno della settimana 
format(d,"%A")

format(d,"%b") # mese
format(d,"%B")

format(d,"%c") # data completa

format(d,"%D") # aa/gg/mm
format(d,"%T") # hh:mm:ss

?ISOdate

Conversioni a POSIX*

d <- date()
d
t <- Sys.time()
t

as.POSIXct(d) # errore
as.POSIXct(t)
as.POSIXlt(t)

as.POSIXlt(t,"GMT") # tiene conto del fuso orario
as.POSIXlt(t, "Australia/Darwin") # potrebbe non funzionre su alcuni OS

unlist(as.POSIXlt(t)) # scomposizione della data

d
strptime(d,"%A %B  %d %H:%M:%S %Y") # NA

d1 <- format(Sys.time(), "%a %b %d %H:%M:%S %Y") # quasi come date() ma nel locale
d1

strptime(d1,"%a %b %d %H:%M:%S %Y")

Sys.getlocale()
?Sys.setlocale # non su tutti gli OS funziona

months(t)
weekdays(t)
quarters(t)

julian(t) # numero di giorni dal 1 gennaio 1970

.leap.seconds

Conti con le date

ISOdate(2006,6,9) - ISOdate(2006, 3, 1) #100 giorni

as.POSIXct("2006-06-09") - as.POSIXct("2006-03-01") #99.9... e' cambiata l'ora

# infatti CEST, CET la "S" = daylight "S"aving
c(as.POSIXct("2006-06-09"), as.POSIXct("2006-03-01"))

# se l'ora non viene specificata allora
# ISOdate ha come riferimento mezzogiorno
# POSIXct la mezzanotte, 
format(ISOdate(2006,6,9),"%H:%M:%S")
format(as.POSIXct("2006-06-09"),"%H:%M:%S")

# allora, se non ci serve una precisione in secondi, arrotondiamo

round(as.POSIXct("2006-06-09") - as.POSIXct("2006-03-01")) #100 giorni

# poiche' tutte queste date sono memorizzate come numero di secondi dal
# 1 gennaio 1970 o qualche altra data di riferimento

Oggetti ts

split.screen(c(2,2))
t1 <- ts(rnorm(10), start=c(2006,1,1), frequency=1) # 1 anno a partire dal 1 gennaio 2006
t1
screen(1) # scegliamo in qualse screen disegnare
plot(t1)

t2 <- ts(rnorm(10), start=c(2006,1,1), frequency=7) # 7 giorni
t2
screen(2)
plot(t2)

t3 <- ts(rnorm(10), start=c(2006,1,1), frequency=4) # 4 quarti
t3
screen(3)
plot(t3)

t4 <- ts(rnorm(10), start=c(2006,1,1), frequency=12) # 12 mesi
t4
screen(4)
plot(t4)

close.screen(all=TRUE) # rimettiamo a posto il device

# oppure
split.screen(c(2,2))
screen(2)
split.screen(c(1,2))
screen(1); plot(t1)
screen(5); plot(t2)
screen(6); plot(t3)
screen(3); plot(t4)
close.screen(all=TRUE) # rimettiamo a posto il device

t <- ts(rnorm(18), start=c(2006,5,1), frequency=12) 
t
t[4] <- NA
t

write.table(data.frame(PIL=rnorm(10)), file="serie.txt")
dati <- read.table(file="serie.txt", header=TRUE)
dati
PIL <- ts(dati$PIL, start=c(2006,6,1), frequency=1)
plot(PIL)

require(tseries)
PIL <- read.ts(file="serie.txt", header=TRUE, start=c(2006,6,1), frequency=1)
PIL
plot(PIL)

Funzioni accessorie

start(PIL)
end(PIL)
frequency(PIL)
deltat(PIL)
time(PIL)

window(PIL, 2013, 2017) # estrae una finestra di dati

window(PIL, 2013, 2027)
window(PIL, 2013, 2027,extend=T) # aggiungiamo termini alla serie

Z ordered observations

require(zoo)
z1.index <- ISOdatetime(2006, rep(1:4,5), sample(28,10), 0,0,0)
z1.index
z1.data <- rnorm(10)
z1 <- zoo(z1.data, z1.index)
plot(z1)

z2.index <- as.POSIXct(paste(2006, rep(1:4,5), sample(28,10),sep="-"))
z2.index
z2.data <- sin(2*1:10/pi)
z2 <- zoo(z2.data, z2.index)
plot(z2)

Z.index <- as.Date(sample(12450:12500,10))
Z.data <- matrix(rnorm(30), ncol=3)
colnames(Z.data) <- c("A","B","C")
Z <- zoo(Z.data, Z.index)

Z
summary(Z)

matplot(Z)
matplot(Z,type="b",axes=F,pch=colnames(Z))
axis(1,1:10,labels=as.character(time(Z)))
axis(2)
box()

plot(Z)
plot(Z, plot.type="single", col=2:4)

zr1 <- zooreg(sin(1:9), start=2006, frequency=4)
zr2 <- zoo(sin(1:9), seq(2006,2008, by=1/4), frequency=4)

# merge di serie storiche
rbind(z1,z2)
cbind(z1,z2)

merge(z1,z2)
merge(z1,z2, all=FALSE)

merge(z1,1:length(z1),2)

# operazioni algebriche e confronti
z1+z2
z1 < z2
cumsum(Z)

# attenzione a quello che segue !
idx <- rep(z1.index[1],3)
idx
zoo(1:3,idx)

# funzioni accessorie
index(z1)
z1
tmp <- coredata(z1) 
tmp
coredata(z1) <- 1:length(z1)
z1
coredata(z1) <- tmp
z1

start(Z)
end(Z)
window(Z, start=as.Date("2004-03-01"))
window(Z, end=as.Date("2004-03-01"))
 
# missing values 
z1[sample(1:10,3)] <- NA
z1

na.omit(z1)
na.contiguous(z1)
na.approx(z1)

Analisi serie storiche

# medie mobili sum(X_i, i in -a:a)/(2a+1)

t <- ts(rnorm(30), start=c(2006,1,1), frequency=12) 
t.fil <- filter(t, filter=rep(1/2,2))
plot(t)
lines(t.fil, col="red")
t.fil2 <- filter(t, filter=rep(1/9,9))
lines(t.fil2, col="blue")

t
diff(t)
plot(t)
lines(diff(t), col="red")

diff(t, lag=3)
t[4:30] - t[1:27]

diff(diff(t))
diff(t, difference=2)


# require(ast) # http://sirio.stat.unipd.it/index.php?id=libast
# sfilter(t)

deka <- decompose(co2)
plot(deka)
(m <- HoltWinters(co2))
plot(m)
plot(fitted(m))


plot(m)
points(predict(m, ahead=3),col="blue")

(s <- HoltWinters(co2,seas="mult"))
plot(s)
plot(fitted(s))


plot(s)
points(predict(s, ahead=3),col="blue")

plot(decompose(co2)) # additivo/molt
plot(stl(co2,"per")) # anche periodico


shapiro.test(deka$random) # test di normalita' dei residui

require(tseries)
jarque.bera.test(deka$random) # test di normalita'

# require(lmtest)
# bptest() # test Breusch-Pagan di omoschedasticita'


# correlogramma
acf(co2)
acf(na.omit(deka$random)) # decisamente non buoni!!!


# test di autocorrelazione dei residui
# LB =n(n+2) sum(t in 1:k, r^2(t)/(n-t)
Box.test(deka$random,type="Ljung-Box")

# BP =n sum(t in 1:k, r^2(t)
Box.test(deka$random,type="Box-Pierce")

# oppure test di durbin-watson in lmtest