library(readr)
fatalities <- read_csv("fatalities.csv", col_names = FALSE)
## Parsed with column specification:
## cols(
## X1 = col_integer()
## )
#barplot(fatalities)
fatalities <- ts(fatalities, start = c(2002,1), end=c(2017,12),frequency=12)
#par(mfrow=c(1,2))
plot.ts(fatalities, ylab="Fatalaties", col="orange", lwd=2,
main="Death Count due to Insurgency in Jammu & Kashmir")
plot(decompose(fatalities))
log_fatalities <- log(fatalities)
par(mfrow=c(1,2))
boxplot(fatalities~floor(time(fatalities)), col= "skyblue",
main = "Original Series")
boxplot(log_fatalities~floor(time(log_fatalities)), col= "tomato",
main="Log-transfomred Series")
plot.ts(log_fatalities, ylab="Fatalaties", col="royalblue", lwd=2,
main="Death Count due to Insurgency in Jammu & Kashmir")
monthplot(log_fatalities)
d1log_fatalities <- diff(log_fatalities,1)
plot.ts(d1log_fatalities, ylab="Fatalaties", col="orange", lwd=2,
main="Death Count due to Insurgency in Jammu & Kashmir")
abline(h=0, lty=2, col="darkgreen", lwd=2)
par(mfrow=c(1,2))
acf(diff(d1log_fatalities,12),ylim=c(-1,1),col=c(2,rep(1,11)),lwd=1,lag.max=84)
pacf(diff(d1log_fatalities,12),ylim=c(-1,1),col=c(rep(1,11),2),lwd=1,lag.max=84)
# fitting the ARIMA(1,0,1)
mod <- arima(log_fatalities, order = c(0,1,3), seasonal=list(order=c(2,1,0),period=12))
##### Previsions a llarg termini amb el model complet ######
ultim <- c(2016,12)
pred=predict(mod,n.ahead=12)
pr<-ts(c(tail(log_fatalities,1),pred$pred),start=ultim+c(1,0),freq=12)
se<-ts(c(0,pred$se),start=ultim+c(1,0),freq=12)
#Intervals
tl1<-ts(exp(pr-1.96*se),start=ultim+c(1,0),freq=12)
tu1<-ts(exp(pr+1.96*se),start=ultim+c(1,0),freq=12)
pr1<-ts(exp(pr),start=ultim+c(1,0),freq=12)
ts.plot(fatalities,tl1,tu1,pr1,lty=c(1,2,2,1),col=c(1,4,4,2),xlim=c(ultim[1]-2,ultim[1]+3),type="o",main="Model ARIMA(1,0,1)", ylim = c(0,100))
abline(v=(ultim[1]-2):(ultim[1]+3),lty=3,col=4)
#options(warn = 0)
#install.packages("forecast")
library(forecast)
library(TTR)
insurgency_fit <- HoltWinters(log_fatalities)
plot(insurgency_fit)
insurgency_forecasts <- forecast:::forecast.HoltWinters(insurgency_fit, h=36)
forecast:::plot.forecast(insurgency_forecasts, col = "tomato", lwd=2,
main = "Insurgency Forecast for 2018-20", fcol = "navyblue",
xlab = expression("Machine Intelligence Unit (Project Victor Sierra)"))
print(exp(insurgency_forecasts$mean))
## Jan Feb Mar Apr May Jun Jul
## 2018 17.32707 18.23808 23.22308 24.51403 34.39955 33.94007 34.18904
## 2019 15.47511 16.28875 20.74093 21.89391 30.72283 30.31247 30.53483
## 2020 13.82109 14.54776 18.52409 19.55383 27.43910 27.07259 27.27119
## Aug Sep Oct Nov Dec
## 2018 35.82325 37.69716 30.61578 23.29190 19.27359
## 2019 31.99437 33.66799 27.34349 20.80240 17.21358
## 2020 28.57473 30.06947 24.42094 18.57899 15.37374
- SQRT trnasformation
d1d12serie <- diff(diff(log_fatalities,12),1)
validation=function(model,dades){
s=frequency(get(model$series))
resid=model$residuals
par(mfrow=c(2,2),mar=c(3,3,3,3))
#Residuals plot
plot(resid,main="Residuals")
abline(h=0)
abline(h=c(-3*sd(resid),3*sd(resid)),lty=3,col=4)
#Square Root of absolute values of residuals (Homocedasticity)
scatter.smooth(sqrt(abs(resid)),main="Square Root of Absolute residuals",
lpars=list(col=2))
#Normal plot of residuals
qqnorm(resid)
qqline(resid,col=2,lwd=2)
##Histogram of residuals with normal curve
hist(resid,breaks=20,freq=F)
curve(dnorm(x,mean=mean(resid),sd=sd(resid)),col=2,add=T)
#ACF & PACF of residuals
par(mfrow=c(1,2))
acf(resid,ylim=c(-1,1),lag.max=60,col=c(2,rep(1,s-1)),lwd=1)
pacf(resid,ylim=c(-1,1),lag.max=60,col=c(rep(1,s-1),2),lwd=1)
par(mfrow=c(1,1))
#ACF & PACF of square residuals
par(mfrow=c(1,2))
acf(resid^2,ylim=c(-1,1),lag.max=60,col=c(2,rep(1,s-1)),lwd=1)
pacf(resid^2,ylim=c(-1,1),lag.max=60,col=c(rep(1,s-1),2),lwd=1)
par(mfrow=c(1,1))
#Ljung-Box p-values
par(mar=c(2,2,1,1))
tsdiag(model,gof.lag=7*s)
cat("\n--------------------------------------------------------------------\n")
print(model)
#Stationary and Invertible
cat("\nModul of AR Characteristic polynomial Roots: ",
Mod(polyroot(c(1,-model$model$phi))),"\n")
cat("\nModul of MA Characteristic polynomial Roots: ",
Mod(polyroot(c(1,model$model$theta))),"\n")
#Model expressed as an MA infinity (psi-weights)
psis=ARMAtoMA(ar=model$model$phi,ma=model$model$theta,lag.max=36)
names(psis)=paste("psi",1:36)
cat("\nPsi-weights (MA(inf))\n")
cat("\n--------------------\n")
print(psis[1:20])
#Model expressed as an AR infinity (pi-weights)
pis=-ARMAtoMA(ar=-model$model$theta,ma=-model$model$phi,lag.max=36)
names(pis)=paste("pi",1:36)
cat("\nPi-weights (AR(inf))\n")
cat("\n--------------------\n")
print(pis[1:20])
##Shapiro-Wilks Normality test
print(shapiro.test(resid(mod)))
#Sample ACF vs. Teoric ACF
par(mfrow=c(2,2),mar=c(3,3,3,3))
acf(dades, ylim=c(-1,1) ,lag.max=36,main="Sample ACF")
plot(ARMAacf(model$model$phi,model$model$theta,lag.max=36),ylim=c(-1,1),
type="h",xlab="Lag", ylab="", main="ACF Teoric")
abline(h=0)
#Sample PACF vs. Teoric PACF
pacf(dades, ylim=c(-1,1) ,lag.max=36,main="Sample PACF")
plot(ARMAacf(model$model$phi,model$model$theta,lag.max=36, pacf=T),ylim=c(-1,1),
type="h", xlab="Lag", ylab="", main="PACF Teoric")
abline(h=0)
par(mfrow=c(1,1))
}
validation(mod, d1d12serie)
##
## --------------------------------------------------------------------
##
## Call:
## arima(x = log_fatalities, order = c(0, 1, 3), seasonal = list(order = c(2, 1,
## 0), period = 12))
##
## Coefficients:
## ma1 ma2 ma3 sar1 sar2
## -0.7172 -0.0226 -0.0365 -0.7736 -0.2057
## s.e. 0.0757 0.0952 0.0763 0.0748 0.0723
##
## sigma^2 estimated as 0.2285: log likelihood = -126.03, aic = 264.07
##
## Modul of AR Characteristic polynomial Roots: 1.068106 1.068106 1.068106 1.068106 1.068106 1.068106 1.068106 1.068106 1.068106 1.068106 1.068106 1.068106 1.068106 1.068106 1.068106 1.068106 1.068106 1.068106 1.068106 1.068106 1.068106 1.068106 1.068106 1.068106
##
## Modul of MA Characteristic polynomial Roots: 1.246658 4.68529 4.68529
##
## Psi-weights (MA(inf))
##
## --------------------
## psi 1 psi 2 psi 3 psi 4 psi 5 psi 6
## -0.71719513 -0.02258780 -0.03654094 0.00000000 0.00000000 0.00000000
## psi 7 psi 8 psi 9 psi 10 psi 11 psi 12
## 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 -0.77361069
## psi 13 psi 14 psi 15 psi 16 psi 17 psi 18
## 0.55482982 0.01747416 0.02826846 0.00000000 0.00000000 0.00000000
## psi 19 psi 20
## 0.00000000 0.00000000
##
## Pi-weights (AR(inf))
##
## --------------------
## pi 1 pi 2 pi 3 pi 4 pi 5 pi 6
## -0.71719513 -0.53695665 -0.43784349 -0.35235487 -0.28221802 -0.22636352
## pi 7 pi 8 pi 9 pi 10 pi 11 pi 12
## -0.18159688 -0.14566596 -0.11684433 -0.09372617 -0.07518198 -0.83391751
## pi 13 pi 14 pi 15 pi 16 pi 17 pi 18
## -0.60320462 -0.45419899 -0.36984650 -0.29755312 -0.23835453 -0.19118231
## pi 19 pi 20
## -0.15337180 -0.12302559
##
## Shapiro-Wilk normality test
##
## data: resid(mod)
## W = 0.97226, p-value = 0.0007372
