# Analysis of the Index data #
##############################

# Reading the data (after saving the variavle madad in a csv file)
data <- read.csv(file="Madad.csv")
madad <- ts(data$madad,start=1990,frequency=12)

#####################
# Plotting the data #
#####################
require(forecast)


tsdisplay(madad)
tsdisplay(diff(log(madad)))


#####################
# Fitting the model #
#####################
# Fitting the basic model
fit.0 <- auto.arima(log(madad))
tsdiag(fit.0)
fit.0


# Reducing the model
fit.1 <- arima(log(madad), order=c(4,1,0), 
                seasonal=list(order=c(1,0,1), period=12))
tsdiag(fit.1)
fit.1

fit.2 <- arima(log(madad), order=c(3,1,0), 
                seasonal=list(order=c(1,0,1), period=12))
tsdiag(fit.2)
fit.2

fit.3 <- arima(log(madad), order=c(2,1,0), 
                seasonal=list(order=c(1,0,1), period=12))
tsdiag(fit.3)
fit.3

fit.4 <- arima(log(madad), order=c(2,1,0), 
                seasonal=list(order=c(1,0,0), period=12))
tsdiag(fit.4)
fit.4

fit.5 <- arima(log(madad), order=c(1,1,0), 
                seasonal=list(order=c(1,0,1), period=12))
tsdiag(fit.5)
fit.5

fit.arima <- fit.3


#####################
# Checking forecast #
#####################
# Prediction with our model and past observations
madad.past <- ts(data$madad[-c(length(madad)-0:9)],start=1990,frequency=12)
plot(diff(log(madad)))
points(diff(log(madad.past)),col="red")

# 
fit.past <- arima(log(madad.past), order=c(2,1,0), 
                seasonal=list(order=c(1,0,1), period=12))
plot(forecast(fit.past),include=24,main="log(madad)")
points(log(madad))


#####################
# Forecast with BSM #
#####################
fit.BSM <- StructTS(log(madad.past),type="BSM")
plot(forecast(fit.BSM),include=24,main="log(madad)")
points(log(madad))



###############################
# Fitting a model to trend an #
# season adjusted data        #
###############################
plot(decompose(log(madad)))

plot(stl(log(madad),"per"))
madad.adj <- stl(log(madad),"per")$time.series[,"remainder"]
auto.arima(madad.adj)

fit.1 <- arima(madad.adj, order=c(1,0,1), 
                seasonal=list(order=c(2,0,0), period=12),include.mean=FALSE)
tsdiag(fit.1)
fit.1

fit.2 <- arima(madad.adj, order=c(1,0,1), 
                seasonal=list(order=c(1,0,0), period=12),include.mean=FALSE)
tsdiag(fit.2)
fit.2