library(maps)

library(ggplot2)

#地圖資料

data(world.cities);

head(data1)
mycols <- runif(10,min=1,max=length(colors())) 
summary(data1)
qplot(longitude, latitude, data = data1 ,col=rainbow(1:8)) + borders("world", size = 0.5)
data1<-read.csv("C:\\Users\\Administrator\\desktop\\地理座標.csv")

1.result:

                                                          全球2015-2017恐怖襲擊發生區域熱點圖

2.恐怖襲擊死亡的總人數、凶手死亡人數、受傷總數、凶手受傷人數、抓獲的凶手數量預測分析

       我們經過時間序列分析方法去發現未來的反恐態勢，主要圍繞2015年到2017年事件中所有確認死亡的總人數(nkill)、凶手死亡人數(nkillter)、受傷總數（nwound）、凶手受傷人數（nwoundte）、抓獲的凶手數量（ nperpcap）五個方面，按月分別對每個方面進行實踐序列分析，判斷其所屬時間序列模型的哪一種，擬合時間序列函式，預測一年每個方面每個月可能會發生的數量，並提出相關建議。

圖1-1  2015-2017  nkill、nperpcap、nkiller、nwound、nwoundte的時間序列圖

      首先對總死亡人數(nkill)進行時間序列分析，從圖1-1 中總死亡人數的時間序列圖可看出從2015年到2017年總死亡人數的數量呈現一個遞減的趨勢。

圖1-2 nkill的自相關函式圖和偏自相關函式圖

       對總死亡人數的自相關函式ACF和偏自相關函式PACF分析，發現總死亡人數時間序列的自相關函式並沒有快速減為0，存在明顯的拖尾現象。並且其前將近44%的序列呈遞增趨勢，後面部分呈現遞減的趨勢，說明其時間序列不夠平穩，需要進行差分。總死亡人數偏自相關函式都在其偏自相關值範圍內，不存在拖尾現象，採用R軟體函式ndiffs對是否需要差分進行判定，結果支援初步判定。

圖1-3  diff=1的總死亡人數時間序列時序和相關函式圖

        一階差分後的的總死亡人數時間序列的自相關函式和偏自相關函式都在範圍內且圍繞0值分佈，對總死亡人數時間序列進行一階差分使序列變得平穩。但仍然無法判斷是否需要繼續差分，假設需要進行二次差分，得到diff=2的序列、相關函式如圖1-1所示。

圖1-4 diff=2的總死亡人數時間序列時序和相關函式圖

對比diff=1和diff=2的時間序列圖和偏自相關函式可以發現，總死亡人數時間序列二次差分出現了過度差分的情況，所以，一階差分已經滿足平穩總死亡人數時間序列的需求，不需要進行第二次差分。

將總死亡人數時間序列平穩後，我們可以通過acf和pacf圖來判斷模型的階數，如下表。所示還可以通過模型的AIC和BIC值來選取模型

                              表1-1 ARIMA(p,d,q)模型選擇

 通過一階總死亡人數時間序列自相關函式和片自相關函式圖可以我們選取建立ARIMA(2,1,2)模型，採用極大似然估計對模型引數進行估計，最後估計出來的模型為：

   我們通過引數的顯著性檢驗和殘差的正態性和無關性檢驗。 

       引數的顯著性檢驗

       引數的顯著性檢驗是用估計出的係數除以其的標準差得到的商與T統計量5%的臨界值(1.96)比較，商的絕對值大於1.96，則拒絕原假設，認為係數顯著的不為0，否則認為係數不顯著。

   (1)殘差的正態性檢驗：

      畫出殘差的QQ圖即可判斷殘差的正態性，QQ圖中殘差基本完全落在45°線上即為符合正態性假設。否則模型可能出現錯誤。

殘差的無關性檢驗：也稱為殘差的白噪聲檢驗，由白噪聲的定義可知，殘差應為不相關的序列。常用LB統計量來檢驗殘差，LB統計量服從自由度為K-p-q的卡方分佈，公式如下：

Q=n(n+2)∑(ρ/(n−k))

其中n為樣本容量，k為自由度

        擬合函式服從殘差基本完全落在45°線上的要求。通過白噪聲檢驗，結果中P值=0.035<0.05，殘差為非白噪聲序列，說明殘差中還蘊含資訊，模型可以繼續優化，LB值為從畫出的QQ圖和LB檢驗的結果來看，殘差符合正態性假設且不相關，則認為模型擬合數據比較充分，可以用來進行下一步預測。

圖1-6 nkill真實值$擬合值的時間序列                  圖1-7nkill時間序列擬合誤差

由圖可以看出總死亡人數時間序列擬合曲線和實際曲線重合得很好，擬合曲線與實際曲線相比僅僅是向右平移了一點。也就是說，下一時刻的狀態很大程度上是由該時刻決定的。

由於其他四個指標的分析方法與總死亡人數的分析過程類似，在此就不再重複敘述。

1.2015到2017年全球每月被抓捕的恐怖襲擊者人數(nperpcap)擬合函式分析

全球被抓捕的恐怖襲擊者人數時間序列是具有一定下降趨勢且沒有季節性變動的時間序列，可以判斷其屬於相加模型，我們採用Holt-Winter平滑法來對其進行擬合。

Holt-Winters平滑方法的預測模型為：

模型中at、bt由以下公式決定：

兩個係數和β在0-1範圍之間。預測結果如表1-1所示。

表1-2全球被抓捕的恐怖襲擊者人數的擬合函式的自變數係數和相關引數

其中平滑化依靠當前時間點上的水平、趨勢部分的斜率、季節性部分三個引數來控制，分別表示為 alpha，beta 和 gamma。可以看見序列的水平度為0.13，下降趨勢為0.15。全球被抓捕的恐怖襲擊者人數的擬合函式為：

 2.2015到2017年全球恐怖襲擊分子每月死亡人數(nkiller)擬合函式分析

nkiller 是一個不穩定、既有季節趨勢又有下降趨勢的序列。我們採用Holt-Winters 指數平滑法 。Holt-Winters 指數平滑法適用於一個增長或降低趨勢並存在季節性可被描述成為相加模型的時間序列。

三次指數平滑法相比二次指數平滑，增加了第三個量來描述季節性。由於存在季節性，我們採用累加式季節性指數平滑法，對應的等式為：

其中pi為週期性的分量，代表週期的長度。x{i+h}為模型預測的等式。R軟甲計算得到結果如表1-1所示

表1-3  2015到2017年全球恐怖襲擊分子每月總死亡人數(nkiller)擬合函式

      則2015到2017年全球恐怖襲擊分子每月死亡人數(nkiller)擬合函式為：

3.2015到2017年全球恐怖襲擊分子每月受傷總人數(nwound)擬合函式

      該時間序列與上一個時間序列型別一樣。得到其公式：

       該序列不平穩，與每月總死亡人數一樣需要進行一階差分，採用相同辦法擬合得到其為2階自迴歸模型加1階移動平均模型。赤池資訊aic=453.81，標準誤差分別為0.1739、0.1728、0.0868。該模型還存在誤差項，誤差項為-0.8709，對應標準誤差為2.1269。

4.2015到2017年全球恐怖襲擊分子每月死亡人數(nbound)擬合函式

      該序列呈現下降的趨勢，與每月被抓捕的恐怖襲擊者人數時間序列都是沒有季節性變動的時間序列,採用和每月被抓捕的恐怖襲擊分子人數的時間序列一樣的辦法。

                               表1-4 全球恐怖襲擊分子每月死亡人數(nbound)

        該序列呈現下降的趨勢，與每月被抓捕的恐怖襲擊者人數時間序列都是沒有季節性變動的時間序列,採用和每月被抓捕的恐怖襲擊分子人數的時間序列一樣的辦法。

             表1-5全球恐怖襲擊分子每月死亡人數(nbound)

      得到恐怖襲擊分子沒有死亡人數的序列的水平度為0.37，下降趨勢為0。

        根據上訴五個時間序列的擬合函式，我們得到了2018年1月到2018年10月的全球恐怖襲擊中總死亡人數、總受傷人數、恐怖分子被捕獲的人數、恐怖分子受傷人數以及恐怖分子死亡人數。預測結果如下：

library(zoo)
library(xts)
library(TTR)
library(forecast)
library(xlsx)

data1<-read.xlsx("C:/Users/Administrator/desktop/傷亡等統計.xlsx","Sheet3")
nperpcap<-as.data.frame(data1[,1])
nkill<-as.data.frame(data1[,2])
nkillter<-as.data.frame(data1[,3])
nwound<-as.data.frame(data1[,4])
nwoundte<-as.data.frame(data1[,5])

nperpcap1<-ts(nperpcap)
nkill1<-ts(nkill)
nkiller1<-ts(nkillter)
nwound1<-ts(nwound)
nwoundte1<-ts(nwoundte)

par(mfrow=c(2,3))
plot(nperpcap1,col="blue",ylab="nperpcap")
plot(nkill1,col="red",ylab="nkill")
plot(nkiller1,col="black",ylab="nkiller")
plot(nwound1,col="orange",ylab="nwound")
plot(nwoundte1,col="cyan",ylab="nwoundte")

1、畫出時序圖和自相關圖
opar<-par(mfrow=c(1,2))
acf(nkiller1,main="nkill acf",lag=36) #畫出自相關圖
pacf(nkiller1,main="nkill pacf",lag=36)

ndiffs(nkiller1)

ndata<-diff(nkiller1,1)
ndiffs(ndata)

#model1<-HoltWinters(nwound1,gamma = FALSE)
model1<-arima(ndata,order=c(2,1,2),method="ML")
e<-forecast(model1)

4、模型檢驗
Box.test(model1$residuals,lag=6,type='Ljung')
A=qqnorm(model1$residuals,plot=FALSE)
qqline(model1$residuals)
qqnorm(model1$residuals,plot=FALSE)
Box.test(model1$residuals,type="Ljung-Box")


5、預測
p<-forecast(model1,4)

plot(model1)
lines(p$pred,col="blue")
exp(p$pred,col="blue")

6
indata=data1[,4]
n <- length(indata)
u <- matrix(0,n,1)
for(i in 1:n){
  u[i] <- rnorm(1)
}
indata_fit <- matrix(0,n,1)
indata_fit[1] <- indata[1]
indata_fit[2] <- indata[2]
for(i in 2:(n-1)){
  indata_fit[i+1] <- 0.9711*indata[i] - 0.2664*u[i]-0.0793*u[i-1]+u[i+1]
}

opar<-par(mfrow=c(1,2))
plot(indata,type = "l",col=c("blue"),xlab="nkill true $ nkill fit")
lines(indata_fit,type = "l",col=c("red"))

error <- (indata_fit - indata)^2
summary(error)
plot(error,type = "l",xlab="nkill error",col="blue")

boxplot(error)