# ********************************************************************
# First day:
#
#    * Brief introduction to R
#    * Some statistical functions and graphic tools
#    * Time series analysis using astronomical (univariate and multivariate) datasets
#    * Spectral analysis: periodogram
#
# ********************************************************************

# Dataset: "HIP.dat" dataset of stars observed with the
#          European Space Agency's Hipparcos satellite during the 1990s.
#          It gives a table with 9 columns and 2719 rows giving Hipparcos
#          stars lying between 40 and 50 parsecs from the Sun.


setwd("datasets")            # set working directory
getwd()                      # current directory

dir()                        # all files
list.files(pattern=".dat")   # files with extension *.dat

hip = read.table("HIP.dat", header=T)    # read file with names (first row)

head(hip)                    # dataset  (first rows)
dim(hip)                     # matrix dimension
names(hip)                   # variable names

hip[2,]                      # row 2
hip[1:5,1:6]                 # rows 1 to 5, columns 1 to 6
head(hip$Vmag)               # first elements of Vmag variable

hip = hip[,-1]               # removinf the first column (without any statistical value)

# *****************************************************
summary(hip)                 # numerical summaries (all variables)

# Central tendency

apply(hip,2,mean)            # means by columns
apply(hip,2,median)          # median by columns

# ****************************************************
# Dispersion
apply(hip,2,min)            # minimum
apply(hip,2,max)            # maximum
apply(hip,2,range)          # range
apply(hip,2,var)            # variance
apply(hip,2,sd)             # standard deviation

# ****************************************************
# selected function
func <- function(x){ c(min(x), max(x),range(x),var(x)) }
apply(hip,2,func)

# ****************************************************
# Missing values: example
x = c(1:5,NA)               # vector with missing value
x
is.na(x)                    # checks for presence of NAs
any(is.na(x))               # TRUE if there is at least one TRUE

# For dataset...
hip[,9]
any(is.na(hip[,9]))

hip1 = na.omit(hip)         # omit missing values
dim(hip)
dim(hip1)

apply(hip,2,mean)
apply(hip1,2,mean)

apply(hip,2,mean,na.rm=T)   # ignore the NAs without omitting them from the dataset
summary(hip,na.rm=T)

# ****************************************************
attach(hip)                 # to access the individual variables

# Graphics
boxplot(hip)                # parallel boxplots (for all variables)
boxplot(pmRA)               # individual

plot(hip,pch=".")           # scatter plot

plot(RA,DE,xlab="Right ascension",ylab="Declination",
     main="RA and DE from Hipparcos data",col="blue")  # scatter plot for two variables
     
hist(B.V)                   # Histograms
hist(B.V,prob=T)            # # Histograms (probability)

old = par(mfrow=c(1,2))     # print 2 chart in same window
hist(B.V)
hist(B.V,prob=T)

par(old)                    # original setup

hist(B.V)                   # first plot
x11()                       # new window
hist(B.V,prob=T)            # second plot

# ********************************************************************
# Extranting a subset of the data: Hyades stars

# Restrictions:
#    * RA in the range (50,100)
#    * DE in the range (0,25)
#    * pmRA in the range (90,130)
#    * pmDE in the range (-60,-10)
#    * e_Plx <5
#    * Vmag >4 OR B.V <0.2 (this eliminates 4 red giants)

attach(hip)

filter1 =  (RA>50 & RA<100 & DE>0 & DE<25 & pmRA>90 & pmRA<130 & pmDE>-60 & pmDE< -10)
filter2  =  filter1 & e_Plx<5
HyadFilter = filter2 & (Vmag>4 |B.V <0.2)

hyades = hip[HyadFilter,]        # Filtered dataset
dim(hyades)                      # dimension

par(mfrow=c(1,2))
plot(hyades$Vmag,hyades$B.V)     # scatterplot of filtered dataset
plot(hip$Vmag,hip$B.V)           # scatterplot: original dataset

write.table(hyades,file="hyades1.dat",row.names = F)
write(t(hyades),file="hyades2.dat",ncol=8)     #

# **********************************************************
# Time series analysis

# Univariate time series
x = cumsum(rnorm(50))                         # data
y = ts(x, frequency = 12, start = c(1959, 7)) # TS object, start date: July 1959
plot(y)                                       # time series plot


# Multivariate
x = matrix(rnorm(300), 100, 3)                # data matrix
dim(x)                                        # dimension 100 x 3
z <- ts(x, start=c(1961, 1), frequency=12)    # TS object, start date: January 1961
plot(z)                                       # plotted separately
x11()
plot(z, plot.type="single", col=1:3,lwd=2)    # single plot

# Example:
dat   = read.table("sspot.dat",head=T)
sspot = ts(dat,freq=12,start=c(1749,1))

plot(sspot, main = "sunspots data", xlab = "Year",
     ylab = "Monthly sunspot numbers")        # time series plot

# periodic pattern: determine the period using Fourier analysis
sp = spectrum(sspot,plot=T)
names(sp)

plot(sp$freq,sp$spec,type="h",xlab="Frequency") # a visual spectrum

# period: the reciprocal of the frequency associated with maximum spectrum
index = which.max(sp$spec)                      # index of the first maximum
1/sp$freq[index]                                # period (every 11 years)

# Related functions
ffourier = fft(sspot)                           # Fast Discrete Fourier Transform
f        = 0:(nss)/nss                          # frequency 1/t
I        = abs(fft(sspot)/sqrt(nss))^2          # Classical periodogram
plot(f[2:2000],I[2:2000],type="h",xlab="Frequency")


par(mfrow=c(2,1))
acf(sspot)                                      # autocorrelation function
pacf(sspot)                                     # partial autocorrelation






