par(mfrow=c(1,2))
plot(vitelli$temperature, type='l', ylab = "temperatura"); grid()
boxplot(vitelli$temperature, main = "temperatura", col = "palegreen4")
Giorgio Garziano
Il dataset “vitelli” e’ stato analizzato allo scopo di verificare relazioni tra le serie temporali basate sulle misure storiche di ammoniaca, acido solfidrico e diossido di carbonio.
Dal punto di vista dell’analisi prendendo ogni serie temporale singolarmente:
Dal punto di vista invece dell’analisi che considera le relazioni tra le serie temporali:
Si puo’ definire un modello di predizione dei picchi positivi dell’ammoniaca in funzione di valori passati di ammoniaca e diossido di carbonio.
I valori passati di ammoniaca giocano un ruolo importante per la previsione dei picchi di ammoniaca stessa.
analisi esplorativa dei dati
individuazione delle caratteristiche di base delle serie temporali
modellizzazione delle serie temporali una ad una (analisi univariata)
modellizzazione delle serie temporali in relazione le une con le altre (analisi multivariata)
definizione dei modelli di “machine learning” per la predizione di picchi di ammoniaca
clean-up dei dati
sommarizzazione e estrazione di statistiche di base
line-plot, box-plot e density-plot
decomposizione additiva delle serie temporali
test per la presenza di effetti eteroschedatici
test per la presenza delle radici unitarie
individuazione delle variazioni di livello (“structural breaks”)
individuazione dei picchi
modellizzazione ARMA-GARCH
calcolo delle cross-correlazioni
modellizzazione VAR
analisi della stabilità
risposta all’impulso
causalità di Granger
test di cointegrazione
modellizzazione VECM
dataset con i dati storici e picchi di ammoniaca
definizione del training e del test dataset
configurazione delle opzioni di training
fit del modello supervised learning basato su regressione logistica
valutazione dell’importanza dei predittori
ROC plot
tuning della soglia di decisione
valutazione della “confusion matrix”
Si e’ fatto uso del linguaggio R in ambiente RStudio su Windows. Le funzionalità utilizzate sono fornite da 16 librerie R (“packages”).
Questa presentazione viene generata per mezzo della stesura di un file R markdown composto sia da codice che testo. I grafici sono generati dal codice incluso.
La compilazione del file R markdown produce un file HTML che puo’ essere visualizzato in modalità di presentazione a slide.
vitelli <- read.csv("vitelli4.csv", header =T, stringsAsFactors = F, sep=",")
head(vitelli)## Timestamp ISO.Date temperature humidity ammonia
## 1 1.51835e+12 2018-02-11T13:03:48.013Z 10.5 63.6 11
## 2 1.51835e+12 2018-02-11T13:03:48.014Z 10.5 63.6 11
## 3 1.51836e+12 2018-02-11T13:18:47.939Z 10.6 62.7 10
## 4 1.51836e+12 2018-02-11T13:18:47.944Z 10.6 62.8 10
## 5 1.51836e+12 2018-02-11T13:33:48.419Z 10.6 63.8 10
## 6 1.51836e+12 2018-02-11T13:33:48.420Z 10.6 63.5 11
## methane hydrogensulfide carbondioxide lux security
## 1 0 0 554 714 1023
## 2 0 0 558 738 1023
## 3 0 0 548 763 1023
## 4 0 0 548 763 1023
## 5 0 0 541 489 1023
## 6 0 0 548 501 1023vitelli <- vitelli[,-1]sel_idx <- compute_duplicate_dates()
head(vitelli[sel_idx,], 20)## ISO.Date temperature humidity ammonia methane
## 1 2018-02-11T13:03:48.013Z 10.5 63.6 11 0
## 3 2018-02-11T13:18:47.939Z 10.6 62.7 10 0
## 5 2018-02-11T13:33:48.419Z 10.6 63.8 10 0
## 7 2018-02-11T13:48:48.357Z 10.5 64.4 11 0
## 9 2018-02-11T14:03:47.871Z 10.6 62.9 11 0
## 10 2018-02-11T14:18:47.294Z 10.6 63.8 11 0
## 12 2018-02-11T14:33:48.456Z 10.5 63.2 10 0
## 14 2018-02-11T14:48:47.930Z 10.7 64.2 11 0
## 16 2018-02-11T15:03:47.928Z 10.6 63.6 11 0
## 18 2018-02-11T15:18:47.654Z 10.6 63.0 11 0
## 19 2018-02-11T15:33:48.518Z 10.6 64.3 11 0
## 21 2018-02-11T15:48:47.932Z 10.6 64.0 12 0
## 23 2018-02-11T16:03:48.580Z 10.6 64.1 11 0
## 25 2018-02-11T16:18:47.751Z 10.6 63.4 11 0
## 27 2018-02-11T16:33:48.190Z 10.5 65.1 11 0
## 28 2018-02-11T16:48:48.734Z 10.6 66.8 12 0
## 30 2018-02-11T17:03:48.272Z 11.0 66.9 12 0
## 32 2018-02-11T17:18:48.120Z 11.0 68.2 13 0
## 34 2018-02-11T17:33:48.095Z 10.9 69.1 13 0
## 36 2018-02-11T17:48:48.639Z 10.7 69.5 15 0
## hydrogensulfide carbondioxide lux security
## 1 0 554 714 1023
## 3 0 548 763 1023
## 5 0 541 489 1023
## 7 0 548 750 1023
## 9 0 545 738 1023
## 10 0 551 692 1023
## 12 0 558 582 1023
## 14 0 561 519 1023
## 16 0 551 467 1023
## 18 0 554 415 1023
## 19 0 561 366 1023
## 21 0 574 289 1023
## 23 0 577 195 1023
## 25 0 561 78 1023
## 27 0 567 29 1023
## 28 0 580 76 1023
## 30 0 622 74 1023
## 32 0 648 2 1023
## 34 0 648 2 1023
## 36 0 670 2 1022tail(vitelli[sel_idx,], 20)## ISO.Date temperature humidity ammonia methane
## 22967 2018-07-28T03:16:27.081Z 23.2 70.5 NA 0
## 22968 2018-07-28T03:31:27.601Z 23.2 70.3 NA 0
## 22969 2018-07-28T03:46:27.181Z 23.0 70.5 NA 0
## 22971 2018-07-28T04:01:27.241Z 23.1 71.0 NA 0
## 22972 2018-07-28T04:16:27.339Z 22.9 70.6 NA 0
## 22974 2018-07-28T04:31:27.961Z 22.9 71.3 NA 0
## 22976 2018-07-28T04:46:27.281Z 23.0 71.7 NA 0
## 22977 2018-07-28T05:01:27.300Z 23.8 70.7 NA 0
## 22979 2018-07-28T05:16:27.342Z 25.8 66.3 NA 0
## 22980 2018-07-28T05:31:27.661Z 25.5 66.7 NA 0
## 22982 2018-07-28T05:46:27.424Z 24.8 70.0 NA 0
## 22983 2018-07-28T06:01:27.400Z 25.0 68.0 NA 0
## 22984 2018-07-28T06:16:27.523Z 25.4 67.2 NA 0
## 22985 2018-07-28T06:31:27.762Z 25.6 65.4 NA 0
## 22986 2018-07-28T06:46:27.521Z 25.9 64.9 NA 0
## 22987 2018-07-28T07:01:27.560Z 26.3 64.8 NA 0
## 22988 2018-07-28T07:16:27.621Z 26.7 64.6 NA 0
## 22990 2018-07-28T07:31:27.804Z 27.1 63.3 NA 0
## 22991 2018-07-28T08:01:27.846Z NA NA NA NA
## 22992 2018-07-28T08:31:28.061Z NA NA NA NA
## hydrogensulfide carbondioxide lux security
## 22967 0 561 4 1023
## 22968 0 561 9 1023
## 22969 0 570 36 1023
## 22971 0 564 111 1023
## 22972 0 567 316 1023
## 22974 0 570 532 1023
## 22976 0 567 881 1023
## 22977 0 564 2083 1023
## 22979 0 570 1602 1023
## 22980 0 551 1602 1023
## 22982 0 541 1716 1023
## 22983 0 532 1846 1023
## 22984 0 535 1998 1023
## 22985 0 532 1998 1023
## 22986 0 532 2083 1023
## 22987 0 545 2176 1023
## 22988 0 538 2277 1023
## 22990 0 551 2277 1023
## 22991 NA NA NA NA
## 22992 NA NA NA NAvitelli <- vitelli[sel_idx,]
dim(vitelli)## [1] 11496 9A seguire test applicati alle singole serie temporali considerandole una ad una.
basicStats(vitelli[,"ammonia"])## X..vitelli..ammonia
## nobs 11496.000000
## NAs 5376.000000
## Minimum 0.000000
## Maximum 55.000000
## 1. Quartile 5.000000
## 3. Quartile 8.000000
## Mean 7.426748
## Median 6.000000
## Sum 45451.700000
## SE Mean 0.052268
## LCL Mean 7.324285
## UCL Mean 7.529212
## Variance 16.719557
## Stdev 4.088956
## Skewness 2.211014
## Kurtosis 8.234634basicStats(vitelli[,"hydrogensulfide"])## X..vitelli..hydrogensulfide
## nobs 11496.000000
## NAs 117.000000
## Minimum 0.000000
## Maximum 441.000000
## 1. Quartile 0.000000
## 3. Quartile 0.000000
## Mean 0.936198
## Median 0.000000
## Sum 10653.000000
## SE Mean 0.131505
## LCL Mean 0.678426
## UCL Mean 1.193971
## Variance 196.783238
## Stdev 14.027945
## Skewness 25.079925
## Kurtosis 716.289731basicStats(vitelli[,"carbondioxide"])## X..vitelli..carbondioxide
## nobs 1.149600e+04
## NAs 1.170000e+02
## Minimum 4.250000e+02
## Maximum 9.730000e+02
## 1. Quartile 5.380000e+02
## 3. Quartile 5.700000e+02
## Mean 5.592795e+02
## Median 5.480000e+02
## Sum 6.364041e+06
## SE Mean 3.068290e-01
## LCL Mean 5.586780e+02
## UCL Mean 5.598809e+02
## Variance 1.071262e+03
## Stdev 3.273014e+01
## Skewness 2.644670e+00
## Kurtosis 1.244439e+01basicStats(vitelli[,"methane"])## X..vitelli..methane
## nobs 11496
## NAs 117
## Minimum 0
## Maximum 0
## 1. Quartile 0
## 3. Quartile 0
## Mean 0
## Median 0
## Sum 0
## SE Mean 0
## LCL Mean 0
## UCL Mean 0
## Variance 0
## Stdev 0
## Skewness NaN
## Kurtosis NaNbasicStats(vitelli[,"temperature"])## X..vitelli..temperature
## nobs 11496.000000
## NAs 117.000000
## Minimum 3.500000
## Maximum 32.800000
## 1. Quartile 11.500000
## 3. Quartile 24.700000
## Mean 18.663749
## Median 20.100000
## Sum 212374.800000
## SE Mean 0.068112
## LCL Mean 18.530237
## UCL Mean 18.797261
## Variance 52.790615
## Stdev 7.265715
## Skewness -0.191441
## Kurtosis -1.226906basicStats(vitelli[,"humidity"])## X..vitelli..humidity
## nobs 11496.000000
## NAs 117.000000
## Minimum 32.800000
## Maximum 92.000000
## 1. Quartile 61.300000
## 3. Quartile 77.100000
## Mean 68.329423
## Median 69.700000
## Sum 777520.500000
## SE Mean 0.104116
## LCL Mean 68.125336
## UCL Mean 68.533509
## Variance 123.350989
## Stdev 11.106349
## Skewness -0.568060
## Kurtosis -0.247551basicStats(vitelli[,"lux"])## X..vitelli..lux
## nobs 1.149600e+04
## NAs 1.170000e+02
## Minimum 0.000000e+00
## Maximum 5.639000e+03
## 1. Quartile 4.000000e+00
## 3. Quartile 1.063000e+03
## Mean 5.671657e+02
## Median 2.180000e+02
## Sum 6.453779e+06
## SE Mean 6.351051e+00
## LCL Mean 5.547166e+02
## UCL Mean 5.796149e+02
## Variance 4.589816e+05
## Stdev 6.774818e+02
## Skewness 1.048104e+00
## Kurtosis 4.269390e-01par(mfrow=c(1,2))
plot(vitelli$temperature, type='l', ylab = "temperatura"); grid()
boxplot(vitelli$temperature, main = "temperatura", col = "palegreen4")
par(mfrow=c(1,2))
plot(vitelli$humidity, type='l', ylab = "umidità"); grid()
boxplot(vitelli$humidity, main = "umidità", col = "palegreen4")
par(mfrow=c(1,2))
plot(vitelli$lux, type='l', ylab = "luminosità"); grid()
boxplot(vitelli$lux, main = "luminosità", col = "palegreen4")
par(mfrow=c(1,2))
plot(vitelli$ammonia, type='l', ylab = "ammoniaca"); grid()
boxplot(vitelli$ammonia, main = "ammonia", col = "palegreen4")
par(mfrow=c(1,2))
plot(vitelli$hydrogensulfide, type='l', ylab = "acido solfidrico"); grid()
boxplot(vitelli$ammonia, main = "acido solfidrico", col = "palegreen4")
par(mfrow=c(1,2))
plot(vitelli$carbondioxide, type='l', ylab = "diossido di carbonio"); grid()
boxplot(vitelli$carbondioxide, main = "diossido di carbonio", col = "palegreen4")
par(mfrow=c(1,2))
plot(vitelli$methane, type='l', ylab = "metano"); grid()
boxplot(vitelli$methane, main = "metano", col = "palegreen4")
var_of_interest <- c("ISO.Date", "ammonia", "hydrogensulfide", "carbondioxide")
vitelli <- vitelli[,var_of_interest]
vitelli <- vitelli[complete.cases(vitelli),]
dim(vitelli)## [1] 6120 4vitelli <- vitelli[-(1:41),]
head(vitelli)## ISO.Date ammonia hydrogensulfide carbondioxide
## 83 2018-02-12T00:03:47.495Z 14 0 677
## 85 2018-02-12T00:18:46.789Z 13 0 677
## 87 2018-02-12T00:33:48.071Z 14 0 674
## 89 2018-02-12T00:48:48.216Z 14 0 683
## 91 2018-02-12T01:03:48.944Z 15 0 699
## 93 2018-02-12T01:18:47.624Z 14 0 686par(mfrow=c(1,3))
boxplot(vitelli[,"ammonia"], main = "ammonia", col = "palegreen4")
boxplot(vitelli[,"hydrogensulfide"], main = "hydrogensulfide", col = "palegreen4")
boxplot(vitelli[,"carbondioxide"], main = "carbondioxide", col = "palegreen4")
t.test(vitelli[,"ammonia"])##
## One Sample t-test
##
## data: vitelli[, "ammonia"]
## t = 141.41, df = 6078, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 7.287375 7.492260
## sample estimates:
## mean of x
## 7.389817t.test(vitelli[,"hydrogensulfide"])##
## One Sample t-test
##
## data: vitelli[, "hydrogensulfide"]
## t = 6.0021, df = 6078, p-value = 2.06e-09
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.9152072 1.8030030
## sample estimates:
## mean of x
## 1.359105t.test(vitelli[,"carbondioxide"])##
## One Sample t-test
##
## data: vitelli[, "carbondioxide"]
## t = 1177, df = 6078, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 559.2848 561.1508
## sample estimates:
## mean of x
## 560.2178par(mfrow=c(1,3))
plot(density(vitelli[,"ammonia"]), main = "ammonia")
plot(density(vitelli[,"hydrogensulfide"]), main = "hydrogensulfide")
plot(density(vitelli[,"carbondioxide"]), main = "carbondioxide")
gas <- vitelli[complete.cases(vitelli), -1]
rownames(gas) <- NULL
head(gas)## ammonia hydrogensulfide carbondioxide
## 1 14 0 677
## 2 13 0 677
## 3 14 0 674
## 4 14 0 683
## 5 15 0 699
## 6 14 0 686plot(gas$ammonia, type='l', main ="ammonia", xlab= "sample", ylab = "ammonia"); grid()
plot(gas$hydrogensulfide, type='l', main = "hydrogensulfide", xlab= "sample", ylab = "hydrogensulfide"); grid()
plot(gas$carbondioxide, type='l', main = "carbondioxide", xlab= "sample", ylab = "carbondioxide"); grid()
gas_n <- data.frame(sample=seq(1:nrow(gas)), ammonia_n = normalize(gas$ammonia, method="range"), hydrogensulfide_n = normalize(gas$hydrogensulfide, method="range"), carbondioxide_n = normalize(gas$carbondioxide, method="range"))
gas_n_wtol <- melt(gas_n, id.vars="sample", measure.vars=c("ammonia_n", "hydrogensulfide_n", "carbondioxide_n"))
ggplot(gas_n_wtol, aes(x=sample, y=value, color = variable)) + geom_line(size = 0.6) + theme_light()
gas_n_wtol <- melt(gas_n[1:1000,], id.vars="sample",
measure.vars=c("ammonia_n", "hydrogensulfide_n", "carbondioxide_n"))
ggplot(gas_n_wtol, aes(x=sample, y=value, color = variable)) + geom_line(size = 0.6) + theme_light()
gas_n_wtol <- melt(gas_n[1001:2000,], id.vars="sample",
measure.vars=c("ammonia_n", "hydrogensulfide_n", "carbondioxide_n"))
ggplot(gas_n_wtol, aes(x=sample, y=value, color = variable)) + geom_line(size = 0.6) + theme_light()
ammonia_ts <- ts(gas$ammonia, frequency = 96)
hydrogensulfide_ts <- ts(gas$hydrogensulfide, frequency = 96)
carbondioxide_ts <- ts(gas$carbondioxide, frequency = 96)ammonia_ts_dec <- decompose(ammonia_ts)
plot(ammonia_ts_dec); grid()
plot(decompose(hydrogensulfide_ts)); grid()
plot(decompose(carbondioxide_ts)); grid()
tsdisplay(ammonia_ts, lag.max = 96); grid()
tsdisplay(hydrogensulfide_ts, lag.max = 96); grid()
tsdisplay(carbondioxide_ts, lag.max = 96); grid()
Scopo dei seguenti test e’ verificare la eteroschedasticità delle serie temporali in esame.
ammonia_mean <- mean(ammonia_ts)
ammonia_ts_demean <- ammonia_ts - ammonia_mean
ArchTest(ammonia_ts_demean)##
## ARCH LM-test; Null hypothesis: no ARCH effects
##
## data: ammonia_ts_demean
## Chi-squared = 2639.7, df = 12, p-value < 2.2e-16hydrogensulfide_mean <- mean(hydrogensulfide_ts)
hydrogensulfide_ts_demean <- hydrogensulfide_ts - hydrogensulfide_mean
ArchTest(hydrogensulfide_ts_demean)##
## ARCH LM-test; Null hypothesis: no ARCH effects
##
## data: hydrogensulfide_ts_demean
## Chi-squared = 633.47, df = 12, p-value < 2.2e-16carbondioxide_mean <- mean(carbondioxide_ts)
carbondioxide_ts_demean <- carbondioxide_ts - carbondioxide_mean
ArchTest(carbondioxide_ts_demean)##
## ARCH LM-test; Null hypothesis: no ARCH effects
##
## data: carbondioxide_ts_demean
## Chi-squared = 4225.3, df = 12, p-value < 2.2e-16Scopo dei seguenti test e’ verificare la l’ordine di integrazione (“unit-root”) e la natura stazionario delle serie temporali in esame.
(urdftest_lag = floor(12* (nrow(gas)/100)^0.25))## [1] 33summary(ur.df(ammonia_ts, type = "none", lags = urdftest_lag, selectlags="BIC"))##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression none
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -18.089 -0.538 -0.025 0.599 36.265
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 -0.005462 0.002590 -2.109 0.035 *
## z.diff.lag1 -0.725104 0.013007 -55.748 < 2e-16 ***
## z.diff.lag2 -0.602198 0.015967 -37.715 < 2e-16 ***
## z.diff.lag3 -0.471886 0.017679 -26.691 < 2e-16 ***
## z.diff.lag4 -0.381059 0.018574 -20.516 < 2e-16 ***
## z.diff.lag5 -0.304865 0.019025 -16.024 < 2e-16 ***
## z.diff.lag6 -0.235777 0.019179 -12.294 < 2e-16 ***
## z.diff.lag7 -0.202059 0.019004 -10.632 < 2e-16 ***
## z.diff.lag8 -0.157424 0.018526 -8.498 < 2e-16 ***
## z.diff.lag9 -0.111469 0.017601 -6.333 2.58e-10 ***
## z.diff.lag10 -0.064748 0.015843 -4.087 4.43e-05 ***
## z.diff.lag11 -0.088672 0.012785 -6.936 4.47e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.665 on 6033 degrees of freedom
## Multiple R-squared: 0.3566, Adjusted R-squared: 0.3554
## F-statistic: 278.7 on 12 and 6033 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -2.1088
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.58 -1.95 -1.62summary(ur.df(ammonia_ts, type = "drift", lags = urdftest_lag, selectlags="BIC"))##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression drift
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.920 -0.574 -0.066 0.556 36.438
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.170156 0.046947 3.624 0.000292 ***
## z.lag.1 -0.023785 0.005679 -4.188 2.85e-05 ***
## z.diff.lag1 -0.709147 0.013719 -51.689 < 2e-16 ***
## z.diff.lag2 -0.587913 0.016431 -35.782 < 2e-16 ***
## z.diff.lag3 -0.459120 0.018009 -25.493 < 2e-16 ***
## z.diff.lag4 -0.369627 0.018822 -19.638 < 2e-16 ***
## z.diff.lag5 -0.294669 0.019213 -15.337 < 2e-16 ***
## z.diff.lag6 -0.226734 0.019321 -11.735 < 2e-16 ***
## z.diff.lag7 -0.194160 0.019110 -10.160 < 2e-16 ***
## z.diff.lag8 -0.150753 0.018599 -8.106 6.30e-16 ***
## z.diff.lag9 -0.106112 0.017646 -6.014 1.92e-09 ***
## z.diff.lag10 -0.060870 0.015863 -3.837 0.000126 ***
## z.diff.lag11 -0.086435 0.012787 -6.760 1.51e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.663 on 6032 degrees of freedom
## Multiple R-squared: 0.358, Adjusted R-squared: 0.3568
## F-statistic: 280.3 on 12 and 6032 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -4.1882 8.7963
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.43 -2.86 -2.57
## phi1 6.43 4.59 3.78summary(ur.df(ammonia_ts, type = "trend", lags = urdftest_lag, selectlags="BIC"))##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression trend
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.871 -0.577 -0.076 0.551 36.491
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.865e-01 9.350e-02 4.134 3.61e-05 ***
## z.lag.1 -3.582e-02 7.242e-03 -4.946 7.78e-07 ***
## tt -4.183e-05 1.563e-05 -2.676 0.007479 **
## z.diff.lag1 -6.984e-01 1.429e-02 -48.890 < 2e-16 ***
## z.diff.lag2 -5.782e-01 1.682e-02 -34.372 < 2e-16 ***
## z.diff.lag3 -4.503e-01 1.830e-02 -24.604 < 2e-16 ***
## z.diff.lag4 -3.616e-01 1.905e-02 -18.982 < 2e-16 ***
## z.diff.lag5 -2.874e-01 1.939e-02 -14.821 < 2e-16 ***
## z.diff.lag6 -2.202e-01 1.946e-02 -11.315 < 2e-16 ***
## z.diff.lag7 -1.884e-01 1.922e-02 -9.804 < 2e-16 ***
## z.diff.lag8 -1.459e-01 1.868e-02 -7.810 6.72e-15 ***
## z.diff.lag9 -1.022e-01 1.770e-02 -5.773 8.17e-09 ***
## z.diff.lag10 -5.802e-02 1.589e-02 -3.651 0.000263 ***
## z.diff.lag11 -8.479e-02 1.280e-02 -6.627 3.73e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.662 on 6031 degrees of freedom
## Multiple R-squared: 0.3588, Adjusted R-squared: 0.3574
## F-statistic: 259.6 on 13 and 6031 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -4.946 8.2565 12.3589
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau3 -3.96 -3.41 -3.12
## phi2 6.09 4.68 4.03
## phi3 8.27 6.25 5.34Phillips-Perron test per la presenza di radici unitarie. Null hypothesis: non stazionarietà.
pp.test(ammonia_ts, alternative="stationary", type="Z(alpha)", lshort = FALSE)##
## Phillips-Perron Unit Root Test
##
## data: ammonia_ts
## Dickey-Fuller Z(alpha) = -3095.4, Truncation lag parameter = 33,
## p-value = 0.01
## alternative hypothesis: stationarypp.test(ammonia_ts, alternative="stationary", type="Z(t_alpha)", lshort=FALSE)##
## Phillips-Perron Unit Root Test
##
## data: ammonia_ts
## Dickey-Fuller Z(t_alpha) = -40.092, Truncation lag parameter = 33,
## p-value = 0.01
## alternative hypothesis: stationaryKwiatkowski-Phillips-Schmidt-Shin (KPSS) test per verificare la stazionarietà della serie temporale. Null Hypothesis: stazionarietà (in livello o trend)
kpss.test(ammonia_ts, null = c("Level"), lshort = FALSE)##
## KPSS Test for Level Stationarity
##
## data: ammonia_ts
## KPSS Level = 4.8969, Truncation lag parameter = 55, p-value = 0.01kpss.test(ammonia_ts, null = c("Trend"), lshort = FALSE)##
## KPSS Test for Trend Stationarity
##
## data: ammonia_ts
## KPSS Trend = 0.71544, Truncation lag parameter = 55, p-value =
## 0.01summary(ur.df(hydrogensulfide_ts, type = "none", lags = urdftest_lag, selectlags="BIC"))##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression none
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -130.1 0.0 0.0 0.0 440.0
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 -0.42833 0.03446 -12.429 < 2e-16 ***
## z.diff.lag1 -0.47621 0.03465 -13.744 < 2e-16 ***
## z.diff.lag2 -0.49583 0.03424 -14.481 < 2e-16 ***
## z.diff.lag3 -0.38223 0.03389 -11.279 < 2e-16 ***
## z.diff.lag4 -0.42909 0.03314 -12.948 < 2e-16 ***
## z.diff.lag5 -0.42031 0.03210 -13.095 < 2e-16 ***
## z.diff.lag6 -0.44174 0.03116 -14.178 < 2e-16 ***
## z.diff.lag7 -0.25664 0.03018 -8.505 < 2e-16 ***
## z.diff.lag8 -0.31849 0.02905 -10.965 < 2e-16 ***
## z.diff.lag9 -0.30359 0.02711 -11.198 < 2e-16 ***
## z.diff.lag10 -0.34236 0.02512 -13.630 < 2e-16 ***
## z.diff.lag11 -0.24143 0.02299 -10.503 < 2e-16 ***
## z.diff.lag12 -0.19678 0.02072 -9.499 < 2e-16 ***
## z.diff.lag13 -0.20118 0.01714 -11.734 < 2e-16 ***
## z.diff.lag14 -0.07103 0.01284 -5.530 3.34e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16.81 on 6030 degrees of freedom
## Multiple R-squared: 0.5016, Adjusted R-squared: 0.5004
## F-statistic: 404.6 on 15 and 6030 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -12.4289
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.58 -1.95 -1.62summary(ur.df(hydrogensulfide_ts, type = "drift", lags = urdftest_lag, selectlags="BIC"))##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression drift
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -129.49 -0.61 -0.61 -0.61 439.39
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.61167 0.22139 2.763 0.00575 **
## z.lag.1 -0.44956 0.03529 -12.739 < 2e-16 ***
## z.diff.lag1 -0.45633 0.03537 -12.902 < 2e-16 ***
## z.diff.lag2 -0.47735 0.03487 -13.690 < 2e-16 ***
## z.diff.lag3 -0.36517 0.03443 -10.607 < 2e-16 ***
## z.diff.lag4 -0.41337 0.03361 -12.300 < 2e-16 ***
## z.diff.lag5 -0.40612 0.03249 -12.500 < 2e-16 ***
## z.diff.lag6 -0.42900 0.03148 -13.628 < 2e-16 ***
## z.diff.lag7 -0.24542 0.03043 -8.064 8.8e-16 ***
## z.diff.lag8 -0.30846 0.02926 -10.543 < 2e-16 ***
## z.diff.lag9 -0.29507 0.02727 -10.820 < 2e-16 ***
## z.diff.lag10 -0.33531 0.02523 -13.289 < 2e-16 ***
## z.diff.lag11 -0.23590 0.02306 -10.229 < 2e-16 ***
## z.diff.lag12 -0.19259 0.02076 -9.277 < 2e-16 ***
## z.diff.lag13 -0.19841 0.01716 -11.559 < 2e-16 ***
## z.diff.lag14 -0.06967 0.01285 -5.423 6.1e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16.8 on 6029 degrees of freedom
## Multiple R-squared: 0.5022, Adjusted R-squared: 0.501
## F-statistic: 405.5 on 15 and 6029 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -12.7389 81.14
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.43 -2.86 -2.57
## phi1 6.43 4.59 3.78##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression trend
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -129.54 -1.11 -0.70 -0.16 438.65
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.4208202 0.4492889 3.162 0.00157 **
## z.lag.1 -0.4613517 0.0357378 -12.909 < 2e-16 ***
## tt -0.0002595 0.0001254 -2.069 0.03854 *
## z.diff.lag1 -0.4452928 0.0357587 -12.453 < 2e-16 ***
## z.diff.lag2 -0.4670995 0.0352092 -13.266 < 2e-16 ***
## z.diff.lag3 -0.3557163 0.0347195 -10.245 < 2e-16 ***
## z.diff.lag4 -0.4046572 0.0338599 -11.951 < 2e-16 ***
## z.diff.lag5 -0.3982484 0.0327019 -12.178 < 2e-16 ***
## z.diff.lag6 -0.4219540 0.0316539 -13.330 < 2e-16 ***
## z.diff.lag7 -0.2392093 0.0305713 -7.825 5.97e-15 ***
## z.diff.lag8 -0.3029212 0.0293714 -10.313 < 2e-16 ***
## z.diff.lag9 -0.2903683 0.0273585 -10.613 < 2e-16 ***
## z.diff.lag10 -0.3314185 0.0252964 -13.101 < 2e-16 ***
## z.diff.lag11 -0.2328534 0.0231028 -10.079 < 2e-16 ***
## z.diff.lag12 -0.1902811 0.0207843 -9.155 < 2e-16 ***
## z.diff.lag13 -0.1968892 0.0171758 -11.463 < 2e-16 ***
## z.diff.lag14 -0.0689181 0.0128485 -5.364 8.45e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16.8 on 6028 degrees of freedom
## Multiple R-squared: 0.5026, Adjusted R-squared: 0.5013
## F-statistic: 380.7 on 16 and 6028 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -12.9093 55.5504 83.3255
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau3 -3.96 -3.41 -3.12
## phi2 6.09 4.68 4.03
## phi3 8.27 6.25 5.34pp.test(hydrogensulfide_ts, alternative="stationary", lshort = FALSE)##
## Phillips-Perron Unit Root Test
##
## data: hydrogensulfide_ts
## Dickey-Fuller Z(alpha) = -8927.5, Truncation lag parameter = 33,
## p-value = 0.01
## alternative hypothesis: stationarypp.test(hydrogensulfide_ts, alternative="stationary", type="Z(t_alpha)", lshort = FALSE)##
## Phillips-Perron Unit Root Test
##
## data: hydrogensulfide_ts
## Dickey-Fuller Z(t_alpha) = -78.661, Truncation lag parameter = 33,
## p-value = 0.01
## alternative hypothesis: stationarykpss.test(hydrogensulfide_ts, null = c("Level"), lshort = FALSE)##
## KPSS Test for Level Stationarity
##
## data: hydrogensulfide_ts
## KPSS Level = 0.63351, Truncation lag parameter = 55, p-value =
## 0.01959kpss.test(hydrogensulfide_ts, null = c("Trend"), lshort = FALSE)##
## KPSS Test for Trend Stationarity
##
## data: hydrogensulfide_ts
## KPSS Trend = 0.088736, Truncation lag parameter = 55, p-value =
## 0.1summary(ur.df(carbondioxide_ts, type = "none", lags = urdftest_lag, selectlags="BIC"))##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression none
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -154.03 -4.55 -0.23 4.35 355.23
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 -0.000206 0.000309 -0.667 0.505
## z.diff.lag1 -0.346686 0.012846 -26.989 < 2e-16 ***
## z.diff.lag2 -0.120066 0.013529 -8.875 < 2e-16 ***
## z.diff.lag3 -0.112944 0.013538 -8.342 < 2e-16 ***
## z.diff.lag4 -0.106912 0.013526 -7.904 3.19e-15 ***
## z.diff.lag5 -0.059750 0.012838 -4.654 3.33e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.47 on 6039 degrees of freedom
## Multiple R-squared: 0.1139, Adjusted R-squared: 0.113
## F-statistic: 129.3 on 6 and 6039 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -0.6665
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.58 -1.95 -1.62summary(ur.df(carbondioxide_ts, type = "drift", lags = urdftest_lag, selectlags="BIC"))##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression drift
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -152.23 -4.82 -1.07 4.31 357.06
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 25.012683 2.811861 8.895 < 2e-16 ***
## z.lag.1 -0.044733 0.005015 -8.920 < 2e-16 ***
## z.diff.lag1 -0.315868 0.013225 -23.884 < 2e-16 ***
## z.diff.lag2 -0.095051 0.013733 -6.921 4.93e-12 ***
## z.diff.lag3 -0.090891 0.013678 -6.645 3.30e-11 ***
## z.diff.lag4 -0.087873 0.013609 -6.457 1.15e-10 ***
## z.diff.lag5 -0.046515 0.012843 -3.622 0.000295 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.39 on 6038 degrees of freedom
## Multiple R-squared: 0.1253, Adjusted R-squared: 0.1245
## F-statistic: 144.2 on 6 and 6038 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -8.9198 39.7892
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.43 -2.86 -2.57
## phi1 6.43 4.59 3.78summary(ur.df(carbondioxide_ts, type = "trend", lags = urdftest_lag, selectlags="BIC"))##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression trend
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -152.81 -4.78 -0.95 4.32 356.56
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 30.5936936 3.2039929 9.549 < 2e-16 ***
## z.lag.1 -0.0525813 0.0054586 -9.633 < 2e-16 ***
## tt -0.0003889 0.0001074 -3.622 0.000295 ***
## z.diff.lag1 -0.3102971 0.0133012 -23.328 < 2e-16 ***
## z.diff.lag2 -0.0904776 0.0137771 -6.567 5.55e-11 ***
## z.diff.lag3 -0.0868462 0.0137101 -6.334 2.55e-10 ***
## z.diff.lag4 -0.0843494 0.0136303 -6.188 6.48e-10 ***
## z.diff.lag5 -0.0440287 0.0128482 -3.427 0.000615 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.37 on 6037 degrees of freedom
## Multiple R-squared: 0.1272, Adjusted R-squared: 0.1262
## F-statistic: 125.7 on 7 and 6037 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -9.6327 30.9513 46.4191
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau3 -3.96 -3.41 -3.12
## phi2 6.09 4.68 4.03
## phi3 8.27 6.25 5.34pp.test(carbondioxide_ts, alternative="stationary", lshort = FALSE)##
## Phillips-Perron Unit Root Test
##
## data: carbondioxide_ts
## Dickey-Fuller Z(alpha) = -654.18, Truncation lag parameter = 33,
## p-value = 0.01
## alternative hypothesis: stationarypp.test(carbondioxide_ts, alternative="stationary", type="Z(t_alpha)", lshort = FALSE)##
## Phillips-Perron Unit Root Test
##
## data: carbondioxide_ts
## Dickey-Fuller Z(t_alpha) = -18.465, Truncation lag parameter = 33,
## p-value = 0.01
## alternative hypothesis: stationarykpss.test(carbondioxide_ts, null = c("Level"), lshort = FALSE)##
## KPSS Test for Level Stationarity
##
## data: carbondioxide_ts
## KPSS Level = 2.6192, Truncation lag parameter = 55, p-value = 0.01kpss.test(carbondioxide_ts, null = c("Trend"), lshort = FALSE)##
## KPSS Test for Trend Stationarity
##
## data: carbondioxide_ts
## KPSS Trend = 0.30682, Truncation lag parameter = 55, p-value =
## 0.01ammonia_bp <- breakpoints(ammonia_ts ~ 1)
summary(ammonia_bp)##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = ammonia_ts ~ 1)
##
## Breakpoints at observation number:
##
## m = 1 1327
## m = 2 1327 3564
## m = 3 1327 2238 3412
## m = 4 1327 2238 3410 4511
## m = 5 1327 2238 3149 4257 5168
##
## Corresponding to breakdates:
##
## m = 1 14(79)
## m = 2 14(79) 38(12)
## m = 3 14(79) 24(30) 36(52)
## m = 4 14(79) 24(30) 36(50) 47(95)
## m = 5 14(79) 24(30) 33(77) 45(33) 54(80)
##
## Fit:
##
## m 0 1 2 3 4 5
## RSS 100898 39886 39319 39130 39040 39084
## BIC 34346 28722 28652 28640 28644 28668plot(ammonia_ts, ylab = "carbondioxide", main = "ammonia level breakpoints"); grid()
lines(fitted(ammonia_bp, breaks = 3), col = "red", lwd=3)
lines(confint(ammonia_bp, breaks = 3))
carbondioxide_ts_bp <- breakpoints(carbondioxide_ts ~ 1)
summary(carbondioxide_ts_bp)##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = carbondioxide_ts ~ 1)
##
## Breakpoints at observation number:
##
## m = 1 995
## m = 2 995 4039
## m = 3 995 2102 3013
## m = 4 995 2102 3013 4205
## m = 5 995 2102 3013 4205 5168
##
## Corresponding to breakdates:
##
## m = 1 11(35)
## m = 2 11(35) 43(7)
## m = 3 11(35) 22(86) 32(37)
## m = 4 11(35) 22(86) 32(37) 44(77)
## m = 5 11(35) 22(86) 32(37) 44(77) 54(80)
##
## Fit:
##
## m 0 1 2 3 4 5
## RSS 8369896 6106561 6010138 5812927 5782437 5781392
## BIC 61205 59306 59227 59041 59027 59043plot(carbondioxide_ts, ylab = "carbondioxide", main = "carbondioxide level breakpoints"); grid()
lines(fitted(carbondioxide_ts_bp, breaks = 4), col = "red", lwd=3)
lines(confint(carbondioxide_ts_bp, breaks = 4))
Identifichiamo la locazione temporale dei picchi per le tre serie temporali in esame.
Allo scopo, definiamo una funzione di utility che ci restituisce la serie temporale priva del trend in essa identificabile.
detrend <- function(a_ts) {
a_ts_dec <- decompose(a_ts)$trend
f <- frequency(a_ts)
beg <- (1:(f/2))
end <- ((length(a_ts) - (f/2)) : length(a_ts))
middle <- ((f/2) + 1) : (length(a_ts) - (f/2) - 1)
ts(c(a_ts[beg], a_ts[middle] - a_ts_dec[middle], a_ts[end]), frequency=f)
}par(mfrow=c(1,2))
plot(density(ammonia_ts)); grid()
plot(density(detrend(ammonia_ts))); grid()
left_thresh <- 0.25; right_thresh <- 0.75
ammonia_outlier_absolute <- getOutliersI(as.vector(ammonia_ts), FLim = c(left_thresh, right_thresh),
distribution = "lognormal")
outlierPlot(gas$ammonia, ammonia_outlier_absolute, mode="qq")
plot(gas$ammonia, type='l', ylab = "ammonia", main = "Ammonia absolute peaks"); grid()
points(x=ammonia_outlier_absolute$iRight, y= gas$ammonia[ammonia_outlier_absolute$iRight], col ='red')
points(x=ammonia_outlier_absolute$iLeft, y=gas$ammonia[ammonia_outlier_absolute$iLeft], col ='blue')
length(ammonia_outlier_absolute$iRight)## [1] 369summary(ammonia_ts[ammonia_outlier_absolute$iRight])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 16.00 18.00 19.00 19.62 20.00 55.00boxplot(ammonia_ts[ammonia_outlier_absolute$iRight], col="palegreen4")
left_thresh <- 0.25; right_thresh <- 0.75
ammonia_outlier_relative <- getOutliersI(as.vector(detrend(ammonia_ts)), rho=c(40,40), FLim = c(left_thresh, right_thresh), distribution = "normal")
outlierPlot(gas$ammonia, ammonia_outlier_relative, mode="qq")
f_2 <- frequency(ammonia_ts)/2; l <- length(ammonia_ts); l_2 <- l - f_2; beg <- which(ammonia_outlier_relative$iRight <= f_2); end <- which(ammonia_outlier_relative$iRight >= l_2)
ammonia_outlier_relative$iRight <- ammonia_outlier_relative$iRight[-c(beg, end)]
plot(gas$ammonia, type='l', ylab = "ammonia", main = "Ammonia relative peaks"); points(x=ammonia_outlier_relative$iRight, y= gas$ammonia[ammonia_outlier_relative$iRight], col ='red'); grid(); points(x=ammonia_outlier_relative$iLeft, y=gas$ammonia[ammonia_outlier_relative$iLeft], col ='blue')
length(ammonia_outlier_relative$iRight)## [1] 163summary(ammonia_ts[ammonia_outlier_relative$iRight])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 7.00 10.00 18.00 16.39 20.00 55.00boxplot(ammonia_ts[ammonia_outlier_relative$iRight], col ="palegreen4")
summary(na.omit(diff(ammonia_outlier_relative$iRight)))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.00 1.00 2.00 36.62 15.00 842.00plot(density(na.omit(diff(ammonia_outlier_relative$iRight))), main = "Relative peaks interarrival time density"); grid()
left_thresh <- 0.25; right_thresh <- 0.75
hydrogensulfide_outlier <- getOutliersI(as.vector(hydrogensulfide_ts), FLim = c(left_thresh, right_thresh),
distribution = "normal")
outlierPlot(gas$hydrogensulfide, hydrogensulfide_outlier, mode="qq")
plot(gas$hydrogensulfide, type='l', ylab = "hydrogensulfide", main = "hydrogensulfide"); grid()
points(x=hydrogensulfide_outlier$iRight, y= gas$hydrogensulfide[hydrogensulfide_outlier$iRight], col ='red')
points(x=hydrogensulfide_outlier$iLeft, y=gas$hydrogensulfide[hydrogensulfide_outlier$iLeft], col ='blue')
left_thresh <- 0.25 ; right_thresh <- 0.75
hydrogensulfide_outlier <- getOutliersI(as.vector(detrend(hydrogensulfide_ts)), FLim = c(left_thresh, right_thresh), distribution = "normal")
outlierPlot(gas$hydrogensulfide, hydrogensulfide_outlier, mode="qq")
f_2 <- frequency(hydrogensulfide_ts)/2; l <- length(hydrogensulfide_ts); l_2 <- l - f_2; beg <- which(hydrogensulfide_outlier$iRight <= f_2); end <- which(hydrogensulfide_outlier$iRight >= l_2)
hydrogensulfide_outlier$iRight <- hydrogensulfide_outlier$iRight[-c(beg, end)]
plot(gas$hydrogensulfide, type='l', ylab = "hydrogensulfide", main = "hydrogensulfide"); grid(); points(x=hydrogensulfide_outlier$iRight, y= gas$hydrogensulfide[hydrogensulfide_outlier$iRight], col ='red'); points(x=hydrogensulfide_outlier$iLeft, y=gas$hydrogensulfide[hydrogensulfide_outlier$iLeft], col ='blue')
left_thresh <- 0.25; right_thresh <- 0.75
carbondioxide_outlier <- getOutliersI(as.vector(carbondioxide_ts), FLim = c(left_thresh, right_thresh), distribution = "lognormal")
outlierPlot(gas$carbondioxide, carbondioxide_outlier, mode="qq")
plot(gas$carbondioxide, type='l', ylab = "carbondioxide", main = "carbondioxide"); grid()
points(x=carbondioxide_outlier$iRight, y= gas$carbondioxide[carbondioxide_outlier$iRight], col ='red')
points(x=carbondioxide_outlier$iLeft, y=gas$carbondioxide[carbondioxide_outlier$iLeft], col ='blue')
left_thresh <- 0.25; right_thresh <- 0.75
carbondioxide_outlier <- getOutliersI(as.vector(detrend(carbondioxide_ts)), FLim = c(left_thresh, right_thresh), distribution = "normal")
outlierPlot(gas$carbondioxide, carbondioxide_outlier, mode="qq")
f_2 <- frequency(carbondioxide_ts)/2; l <- length(carbondioxide_ts); l_2 <- l - f_2; beg <- which(carbondioxide_outlier$iRight <= f_2); end <- which(carbondioxide_outlier$iRight >= l_2)
carbondioxide_outlier$iRight <- carbondioxide_outlier$iRight[-c(beg, end)]
plot(gas$carbondioxide, type='l', ylab = "carbondioxide", main = "carbondioxide"); grid(); points(x=carbondioxide_outlier$iRight, y= gas$carbondioxide[carbondioxide_outlier$iRight], col ='red'); points(x=carbondioxide_outlier$iLeft, y=gas$carbondioxide[carbondioxide_outlier$iLeft], col ='blue')
Vediamo prima in che cosa consiste un modello ARMA ed un modello ARMA-GARCH.
\[ \begin{equation} \begin{cases} y_{t}\ -\ \phi_{1}y_{t-1}\ -\ \phi_{2}y_{t-2}\ +\ ...\ + \phi_{p}y_{t}\ =\ \phi_{0}\ +\ u_{t}\ +\ \theta_{1}u_{t-1}\ +\ \theta_{2}u_{t-2}\ +\ + \theta_{q}u_{t-q} \\ \\ u_{t}\ =\ \sigma\epsilon_{t}\\ \\ \epsilon_{t}\ =\ N(0,1) \end{cases} \end{equation} \]
\[ \begin{equation} \begin{cases} y_{t}\ -\ \phi_{1}y_{t-1}\ -\ \phi_{2}y_{t-2}\ +\ ...\ + \phi_{p}y_{t}\ =\ \phi_{0}\ +\ u_{t}\ +\ \theta_{1}u_{t-1}\ +\ \theta_{2}u_{t-2}\ +\ ...\ + \theta_{q}u_{t-q} \\ \\ u_{t}\ =\ \sigma_{t}\epsilon_{t},\ \ \ \ \ \epsilon_{t}=N(0,1) \\ \\ \sigma_{t}^2\ =\ \omega\ + \alpha_{1}u_{t-1}^2\ +\ \alpha_{2}u_{t-2}^2\ +\ \alpha_{t-q}u_{t-q}^2\ +\ \beta_{1}\sigma_{t-1}^2\ +\ \beta_{2}\sigma_{t-2}^2\ +\ \beta_{p}\sigma_{t-p}^2 \end{cases} \end{equation} \]
xspec = ugarchspec(mean.model = list(armaOrder = c(2, 2), include.mean = TRUE),
variance.model = list(garchOrder = c(1,2), model = 'iGARCH', submodel = NULL),
distribution.model = 'sstd')
ammonia_garch <- ugarchfit(xspec, ammonia_ts, solver="hybrid", fit.control=list(stationarity=1), out.sample=48)
ammonia_garch##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : iGARCH(1,2)
## Mean Model : ARFIMA(2,0,2)
## Distribution : sstd
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 14.026039 0.475180 29.5173 0
## ar1 1.439741 0.000361 3992.2872 0
## ar2 -0.440535 0.000404 -1091.4825 0
## ma1 -0.992664 0.014315 -69.3441 0
## ma2 0.205087 0.012726 16.1153 0
## omega 0.503734 0.067935 7.4150 0
## alpha1 0.653945 0.042505 15.3853 0
## beta1 0.283477 0.054697 5.1827 0
## beta2 0.062578 NA NA NA
## skew 1.113630 0.016817 66.2189 0
## shape 2.702282 0.058952 45.8388 0
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 14.026039 0.419344 33.4476 0e+00
## ar1 1.439741 0.000207 6945.2706 0e+00
## ar2 -0.440535 0.000296 -1488.2384 0e+00
## ma1 -0.992664 0.018511 -53.6257 0e+00
## ma2 0.205087 0.014721 13.9313 0e+00
## omega 0.503734 0.103135 4.8842 1e-06
## alpha1 0.653945 0.062753 10.4209 0e+00
## beta1 0.283477 0.046691 6.0714 0e+00
## beta2 0.062578 NA NA NA
## skew 1.113630 0.016986 65.5623 0e+00
## shape 2.702282 0.065184 41.4563 0e+00
##
## LogLikelihood : -8226.012
##
## Information Criteria
## ------------------------------------
##
## Akaike 2.7312
## Bayes 2.7423
## Shibata 2.7312
## Hannan-Quinn 2.7351
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 3.943 4.707e-02
## Lag[2*(p+q)+(p+q)-1][11] 8.797 2.072e-05
## Lag[4*(p+q)+(p+q)-1][19] 14.046 5.618e-02
## d.o.f=4
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.1391 0.7092
## Lag[2*(p+q)+(p+q)-1][8] 0.2866 0.9993
## Lag[4*(p+q)+(p+q)-1][14] 0.9397 0.9997
## d.o.f=3
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[4] 0.1182 0.500 2.000 0.7310
## ARCH Lag[6] 0.1257 1.461 1.711 0.9837
## ARCH Lag[8] 0.1703 2.368 1.583 0.9984
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 9.37
## Individual Statistics:
## mu 0.03353
## ar1 0.06100
## ar2 0.05575
## ma1 2.43148
## ma2 0.58512
## omega 2.95241
## alpha1 1.43341
## beta1 0.14970
## skew 0.65345
## shape 3.65666
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 2.29 2.54 3.05
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.8259 0.4089
## Negative Sign Bias 0.5601 0.5754
## Positive Sign Bias 1.0352 0.3006
## Joint Effect 1.8894 0.5957
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 271.4 1.411e-46
## 2 30 367.9 2.237e-60
## 3 40 377.7 4.917e-57
## 4 50 586.6 3.406e-93
##
##
## Elapsed time : 2.469538\[ \begin{equation} \begin{cases} a_{t}\ =\ 14.02 + y_{t} \\ \\ y_{t}\ -\ 1.43y_{t-1} +\ 0.44y_{t-2}\ = \ u_{t}\ -\ 0.99u_{t-1} + 0.20u_{t-2} \\ \\ u_{t}\ =\ \sigma_{t}\epsilon_{t},\ \ \ \ \ \epsilon_{t}=N(0,1) \\ \\ \sigma_{t}^2\ =\ 0.50\ + 0.65u_{t-1}^2\ +\ 0.28\sigma_{t-1}^2 +\ 0.06\sigma_{t-2}^2 \end{cases} \end{equation} \]
fit_ts <- ts(fitted(ammonia_garch), frequency = frequency(ammonia_ts))
plot(ammonia_ts, ylab = 'ammonia', lwd = 1, lty ='dashed'); grid(); lines(fit_ts, col = 'red', lty = 'dashed')
legend("topright", c("ammonia", "fitted"), fill=c("black", "red")); mape <- accuracy(ammonia_ts, fit_ts)[,"RMSE"]
mes <- paste("RMSE = ", round(mape, 2), sep=''); text(60, 30, mes)
boxplot(ammonia_ts, fit_ts, main = "Ammonia: actual vs. fitted", col = "palegreen4")
spec = getspec(ammonia_garch);
setfixed(spec) <- as.list(coef(ammonia_garch));
l <- length(ammonia_ts)
forecast = ugarchforecast(ammonia_garch, n.ahead = 24, n.roll=10)
plot(forecast, which=2)
ammonia_sigma <- ammonia_garch@fit$sigma
plot(ammonia_sigma, type='l', "ylab" = "conditional SD", main = "Ammonia Conditional SD"); grid()
xspec = ugarchspec(mean.model = list(armaOrder = c(1,1), include.mean = TRUE),
variance.model = list(garchOrder = c(1,1), model = 'sGARCH',
submodel = NULL),
distribution.model = 'norm')
hydrogensulfide_garch <- ugarchfit(xspec, hydrogensulfide_ts, solver="hybrid")
hydrogensulfide_garch##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,1)
## Mean Model : ARFIMA(1,0,1)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.362161 0.000034 10616.7 0
## ar1 0.047170 0.000009 5495.9 0
## ma1 0.022999 0.000004 5496.0 0
## omega 0.490997 0.000090 5451.7 0
## alpha1 0.050180 0.000014 3709.6 0
## beta1 0.912996 0.000159 5730.4 0
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.362161 NaN NaN NaN
## ar1 0.047170 NaN NaN NaN
## ma1 0.022999 NaN NaN NaN
## omega 0.490997 NaN NaN NaN
## alpha1 0.050180 NaN NaN NaN
## beta1 0.912996 NaN NaN NaN
##
## LogLikelihood : -18884.87
##
## Information Criteria
## ------------------------------------
##
## Akaike 6.2151
## Bayes 6.2218
## Shibata 6.2151
## Hannan-Quinn 6.2174
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.6102 0.4347
## Lag[2*(p+q)+(p+q)-1][5] 1.8637 0.9783
## Lag[4*(p+q)+(p+q)-1][9] 2.3565 0.9637
## d.o.f=2
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.004646 0.9457
## Lag[2*(p+q)+(p+q)-1][5] 0.017298 0.9999
## Lag[4*(p+q)+(p+q)-1][9] 0.026428 1.0000
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.01163 0.500 2.000 0.9141
## ARCH Lag[5] 0.01631 1.440 1.667 0.9990
## ARCH Lag[7] 0.02059 2.315 1.543 1.0000
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: NaN
## Individual Statistics:
## mu NaN
## ar1 NaN
## ma1 NaN
## omega NaN
## alpha1 NaN
## beta1 NaN
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.49 1.68 2.12
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.270986 0.7864
## Negative Sign Bias 0.164492 0.8693
## Positive Sign Bias 0.004407 0.9965
## Joint Effect 0.120793 0.9892
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 54148 0
## 2 30 101378 0
## 3 40 97523 0
## 4 50 149598 0
##
##
## Elapsed time : 1.419778\[ \begin{equation} \begin{cases} h(t)\ = 0.36 + y_{t} \\ \\ y_{t}\ -\ 0.05y_{t-1}\ =\ 0.02u_{t-1}\ +\ u_{t} \\ \\ u_{t}\ =\ \sigma_{t}\epsilon_{t},\ \ \ \ \ \epsilon_{t}=N(0,1) \\ \\ \sigma_{t}^2\ =\ 0.49\ + 0.05u_{t-1}^2\ +\ 0.91\sigma_{t-1}^2 \end{cases} \end{equation} \]
fit_ts <- ts(fitted(hydrogensulfide_garch), frequency = frequency(hydrogensulfide_ts))
plot(hydrogensulfide_ts, lwd=2); grid(); lines(fit_ts, col = 'red', lty='dashed')
legend("topright", c("hydrogensulfide", "fitted"), fill=c("black", "red"))
rmse <- accuracy(hydrogensulfide_ts, fitted(hydrogensulfide_garch))[,"RMSE"];
mes <- paste("RMSE = ", round(rmse, 2), sep=''); text(60, 300, mes)
hydrogensulfide_sigma <- hydrogensulfide_garch@fit$sigma
plot(hydrogensulfide_sigma, type='l', "ylab" = "conditional SD", main = "Hydrogensulfide Conditional SD"); grid()
xspec = ugarchspec(mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
variance.model = list(garchOrder = c(1, 2), model = 'iGARCH',
submodel = "GARCH"),
distribution.model = 'snorm')
carbondioxide_garch <- ugarchfit(xspec, carbondioxide_ts, solver="hybrid", fit.control=list(stationarity=0))
carbondioxide_garch##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : iGARCH(1,2)
## Mean Model : ARFIMA(1,0,1)
## Distribution : snorm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 552.177528 1.596761 345.8110 0
## ar1 0.947687 0.004725 200.5788 0
## ma1 -0.272649 0.016615 -16.4101 0
## omega 0.935525 0.153962 6.0764 0
## alpha1 0.074441 0.005851 12.7229 0
## beta1 0.383380 0.061330 6.2511 0
## beta2 0.542179 NA NA NA
## skew 1.060362 0.012101 87.6259 0
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 552.177528 3.467051 159.2643 0.000000
## ar1 0.947687 0.011439 82.8451 0.000000
## ma1 -0.272649 0.027740 -9.8288 0.000000
## omega 0.935525 0.660770 1.4158 0.156831
## alpha1 0.074441 0.035143 2.1182 0.034157
## beta1 0.383380 0.114611 3.3450 0.000823
## beta2 0.542179 NA NA NA
## skew 1.060362 0.053647 19.7657 0.000000
##
## LogLikelihood : -22782.67
##
## Information Criteria
## ------------------------------------
##
## Akaike 7.4978
## Bayes 7.5056
## Shibata 7.4978
## Hannan-Quinn 7.5005
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 2.278 0.13124
## Lag[2*(p+q)+(p+q)-1][5] 3.290 0.30332
## Lag[4*(p+q)+(p+q)-1][9] 8.473 0.03839
## d.o.f=2
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 5.609 0.01787
## Lag[2*(p+q)+(p+q)-1][8] 6.101 0.22518
## Lag[4*(p+q)+(p+q)-1][14] 6.757 0.53751
## d.o.f=3
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[4] 0.001008 0.500 2.000 0.9747
## ARCH Lag[6] 0.343752 1.461 1.711 0.9342
## ARCH Lag[8] 0.597668 2.368 1.583 0.9734
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 3.0675
## Individual Statistics:
## mu 0.6092
## ar1 0.2858
## ma1 0.9835
## omega 0.3450
## alpha1 0.3245
## beta1 0.2906
## skew 0.2853
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.69 1.9 2.35
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 2.827 4.715e-03 ***
## Negative Sign Bias 1.817 6.923e-02 *
## Positive Sign Bias 2.223 2.622e-02 **
## Joint Effect 22.230 5.844e-05 ***
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 972.2 4.482e-194
## 2 30 986.4 6.439e-189
## 3 40 1036.1 6.313e-192
## 4 50 1067.5 1.626e-191
##
##
## Elapsed time : 0.3960159\[ \begin{equation} \begin{cases} c_{t}\ =\ 552\ +\ y_{t} \\ \\ y_{t}\ -\ 0.95y_{t-1}\ = -\ 0.27u_{t-1}\ +\ u_{t} \\ \\ u_{t}\ =\ \sigma_{t}\epsilon_{t},\ \ \ \ \ \epsilon_{t}=N(0,1) \\ \\ \sigma_{t}^2\ =\ 0.94\ + 0.07u_{t-1}^2\ +\ 0.38\sigma_{t-1}^2\ +\ 0.54\sigma_{t-2}^2 \end{cases} \end{equation} \]
carbondioxide_garch_fit <- ts(fitted(carbondioxide_garch), frequency = frequency(carbondioxide_ts))
plot(carbondioxide_ts, ylab = 'carbondioxide', lwd = 1, lty ='dashed'); grid(); lines(carbondioxide_garch_fit, col = 'red', lty = 'dashed'); legend("topright", c("carbondioxide", "fitted"), fill=c("black", "red"))
mape <- accuracy(carbondioxide_ts, carbondioxide_garch_fit)[,"RMSE"]
mes <- paste("RMSE = ", round(mape, 2), sep=''); text(60, 700, mes)
carbondioxide_sigma <- carbondioxide_garch@fit$sigma
plot(carbondioxide_sigma, type='l', "ylab" = "conditional SD", main = "Carbondioxide Conditional SD"); grid()
Nel seguito, l’analisi ha come oggetto le inter-relazioni tra le serie temporali in esame.
Mostriamo le cross-correlazioni tra le diverse serie temporali in funzione del tempo.
\[ \begin{equation} \begin{cases} \Large \rho_{x, y}(t)\ = \frac{cov(x, y)(t)}{\sigma_{x}(t)\ \sigma_{y}(t)} \\ \\ \Large \sigma_{x}^2 (t+1)\ =\ (1-\theta)x^2(t)\ +\ \theta\ \sigma_{x}^2(t)\\ \\ \Large cov(x,y)(t+1)\ =\ (1-\theta)x(t)y(t)\ +\ \theta\ cov(x,y)(t) \end{cases} \end{equation} \]
Reference: “An Introduction to Analysis of Financial Data with R”, Ruey S. Tsay
covmat <- MTS::EWMAvol(cbind(ammonia_ts, hydrogensulfide_ts))
cov_y1y2 <- covmat$Sigma.t[,2]
cor_y1y2 <- cov_y1y2/sqrt(covmat$Sigma.t[,1]*covmat$Sigma.t[,4])
plot(cor_y1y2, type='l', main="Ammonia & Hydrogensulfide Cross-Correlation", xlab = "time tick"); abline(h=0, col='blue', lty='dotted')
boxplot(cor_y1y2, main = "Ammonia & Hydrogensulfide Cross-Correlation", col = "palegreen4")
summary(cor_y1y2)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.9807 -0.0270 0.2528 0.2432 0.6252 0.9597covmat <- MTS::EWMAvol(cbind(ammonia_ts, carbondioxide_ts))
cov_y1y2 <- covmat$Sigma.t[,2]
cor_y1y2 <- cov_y1y2/sqrt(covmat$Sigma.t[,1]*covmat$Sigma.t[,4])
plot(cor_y1y2, type='l', main="Ammonia & Carbondioxide Cross-Correlation", xlab = "time tick"); abline(h=0, col='blue', lty='dotted')
boxplot(cor_y1y2, main = "Ammonia & Carbondioxide Cross-Correlation", col = "palegreen4")
summary(cor_y1y2)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.7266 0.3023 0.6453 0.5320 0.8348 0.9875covmat <- MTS::EWMAvol(cbind(hydrogensulfide_ts, carbondioxide_ts))
cov_y1y2 <- covmat$Sigma.t[,2]
cor_y1y2 <- cov_y1y2/sqrt(covmat$Sigma.t[,1]*covmat$Sigma.t[,4])
plot(cor_y1y2, type='l', main="Hydrogensulfide & Carbondioxide Cross-Correlation", xlab = "time tick"); abline(h=0, col='blue', lty='dotted')
boxplot(cor_y1y2, main = "Hydrogensulfide & Carbondioxide Cross-Correlation", col = "palegreen4")
summary(cor_y1y2)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.94639 -0.10748 0.08842 0.11885 0.34568 0.97041L’analisi della significatività della cross-correlazione prevede l’uso del “pre-whitening” allo scopo di evitare problemi di cross-correlazione spuria.
prewhiten(x=ammonia_ts, y=hydrogensulfide_ts, ylab="CCF", main ="ammonia & hydrogensulfide")
prewhiten(x=ammonia_ts, y=carbondioxide_ts, ylab="CCF", main ="ammonia & carbondioxide")
prewhiten(x=hydrogensulfide_ts, y=carbondioxide_ts, ylab="CCF", main ="hydrogensulfide & carbondioxide")
\[ \begin{equation} \begin{cases} \large \boldsymbol{y}_{t}\ =\ \boldsymbol{A_{1}}y_{t-1}\ +\ \boldsymbol{A_{2}}y_{t-2}\ +\ ...\ +\ \boldsymbol{A_{p}}y_{t-p}\ +\ \boldsymbol{u_{t}}\\ \\ \\ \large u_{t}\ =\ \sigma\epsilon_{t}\\ \\ \\ \large \epsilon_{t}\ =\ N(0,1)\\ \\ \\ \large \boldsymbol{y}_{t}\ :=\ \begin{bmatrix}ammonia_{t}\\ hydrogensulfide_{t}\\ carbondioxide_{t}\end{bmatrix} \end{cases} \end{equation} \]
VARselect(gas, lag.max = 20, type="both")## $selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 20 20 14 20
##
## $criteria
## 1 2 3 4 5
## AIC(n) 12.31333 12.07611 11.97327 11.92969 11.90188
## HQ(n) 12.31910 12.08534 11.98595 11.94583 11.92148
## SC(n) 12.32994 12.10269 12.00981 11.97620 11.95835
## FPE(n) 222644.08661 175626.34256 158461.49579 151704.20171 147543.12357
## 6 7 8 9 10
## AIC(n) 11.88728 11.84547 11.83840 11.83297 11.82061
## HQ(n) 11.91034 11.87199 11.86838 11.86641 11.85751
## SC(n) 11.95372 11.92188 11.92477 11.92931 11.92691
## FPE(n) 145405.60661 139451.56436 138469.06249 137719.44011 136027.32207
## 11 12 13 14 15
## AIC(n) 11.79529 11.78653 11.78553 11.75809 11.75017
## HQ(n) 11.83565 11.83035 11.83281 11.80883 11.80436
## SC(n) 11.91156 11.91277 11.92173 11.90426 11.90630
## FPE(n) 132626.82316 131470.19590 131338.51746 127783.57183 126775.30228
## 16 17 18 19 20
## AIC(n) 11.74509 11.73347 11.72242 11.72188 11.70887
## HQ(n) 11.80274 11.79458 11.78699 11.78991 11.78035
## SC(n) 11.91119 11.90953 11.90845 11.91788 11.91483
## FPE(n) 126133.09383 124675.52504 123306.00317 123239.26844 121645.61020var_fit <- VAR(gas, p=14, type="both", season = 96)
summary(var_fit)##
## VAR Estimation Results:
## =========================
## Endogenous variables: ammonia, hydrogensulfide, carbondioxide
## Deterministic variables: both
## Sample size: 6065
## Log Likelihood: -61174.488
## Roots of the characteristic polynomial:
## 0.9915 0.9698 0.9212 0.9069 0.9069 0.9029 0.9029 0.8977 0.8977 0.8568 0.8568 0.8477 0.8477 0.8461 0.8461 0.8376 0.8376 0.8235 0.8235 0.8093 0.8093 0.8085 0.7982 0.7982 0.7962 0.7962 0.7916 0.7869 0.7869 0.7832 0.7832 0.7788 0.7788 0.7742 0.7742 0.7731 0.7731 0.7151 0.7151 0.7101 0.4594 0.04134
## Call:
## VAR(y = gas, p = 14, type = "both", season = 96L)
##
##
## Estimation results for equation ammonia:
## ========================================
## ammonia = ammonia.l1 + hydrogensulfide.l1 + carbondioxide.l1 + ammonia.l2 + hydrogensulfide.l2 + carbondioxide.l2 + ammonia.l3 + hydrogensulfide.l3 + carbondioxide.l3 + ammonia.l4 + hydrogensulfide.l4 + carbondioxide.l4 + ammonia.l5 + hydrogensulfide.l5 + carbondioxide.l5 + ammonia.l6 + hydrogensulfide.l6 + carbondioxide.l6 + ammonia.l7 + hydrogensulfide.l7 + carbondioxide.l7 + ammonia.l8 + hydrogensulfide.l8 + carbondioxide.l8 + ammonia.l9 + hydrogensulfide.l9 + carbondioxide.l9 + ammonia.l10 + hydrogensulfide.l10 + carbondioxide.l10 + ammonia.l11 + hydrogensulfide.l11 + carbondioxide.l11 + ammonia.l12 + hydrogensulfide.l12 + carbondioxide.l12 + ammonia.l13 + hydrogensulfide.l13 + carbondioxide.l13 + ammonia.l14 + hydrogensulfide.l14 + carbondioxide.l14 + const + trend + sd1 + sd2 + sd3 + sd4 + sd5 + sd6 + sd7 + sd8 + sd9 + sd10 + sd11 + sd12 + sd13 + sd14 + sd15 + sd16 + sd17 + sd18 + sd19 + sd20 + sd21 + sd22 + sd23 + sd24 + sd25 + sd26 + sd27 + sd28 + sd29 + sd30 + sd31 + sd32 + sd33 + sd34 + sd35 + sd36 + sd37 + sd38 + sd39 + sd40 + sd41 + sd42 + sd43 + sd44 + sd45 + sd46 + sd47 + sd48 + sd49 + sd50 + sd51 + sd52 + sd53 + sd54 + sd55 + sd56 + sd57 + sd58 + sd59 + sd60 + sd61 + sd62 + sd63 + sd64 + sd65 + sd66 + sd67 + sd68 + sd69 + sd70 + sd71 + sd72 + sd73 + sd74 + sd75 + sd76 + sd77 + sd78 + sd79 + sd80 + sd81 + sd82 + sd83 + sd84 + sd85 + sd86 + sd87 + sd88 + sd89 + sd90 + sd91 + sd92 + sd93 + sd94 + sd95
##
## Estimate Std. Error t value Pr(>|t|)
## ammonia.l1 2.310e-01 1.313e-02 17.585 < 2e-16 ***
## hydrogensulfide.l1 3.895e-03 1.268e-03 3.072 0.002133 **
## carbondioxide.l1 9.289e-03 1.609e-03 5.774 8.14e-09 ***
## ammonia.l2 1.255e-01 1.352e-02 9.288 < 2e-16 ***
## hydrogensulfide.l2 -5.331e-03 1.273e-03 -4.187 2.87e-05 ***
## carbondioxide.l2 2.490e-03 1.908e-03 1.306 0.191762
## ammonia.l3 1.468e-01 1.362e-02 10.778 < 2e-16 ***
## hydrogensulfide.l3 1.101e-03 1.273e-03 0.865 0.387132
## carbondioxide.l3 -3.762e-03 1.939e-03 -1.940 0.052378 .
## ammonia.l4 9.891e-02 1.367e-02 7.238 5.12e-13 ***
## hydrogensulfide.l4 -3.399e-04 1.274e-03 -0.267 0.789629
## carbondioxide.l4 -1.059e-03 1.940e-03 -0.546 0.585237
## ammonia.l5 7.828e-02 1.369e-02 5.720 1.12e-08 ***
## hydrogensulfide.l5 -6.602e-04 1.276e-03 -0.518 0.604801
## carbondioxide.l5 -3.580e-03 1.938e-03 -1.847 0.064784 .
## ammonia.l6 5.795e-02 1.371e-02 4.228 2.40e-05 ***
## hydrogensulfide.l6 1.709e-03 1.276e-03 1.340 0.180445
## carbondioxide.l6 -5.533e-04 1.940e-03 -0.285 0.775457
## ammonia.l7 2.426e-02 1.371e-02 1.769 0.076885 .
## hydrogensulfide.l7 2.345e-04 1.274e-03 0.184 0.854019
## carbondioxide.l7 -9.865e-04 1.941e-03 -0.508 0.611229
## ammonia.l8 3.080e-02 1.372e-02 2.245 0.024799 *
## hydrogensulfide.l8 1.136e-03 1.274e-03 0.892 0.372622
## carbondioxide.l8 3.360e-03 1.940e-03 1.732 0.083374 .
## ammonia.l9 4.881e-02 1.370e-02 3.562 0.000372 ***
## hydrogensulfide.l9 7.095e-04 1.275e-03 0.556 0.577994
## carbondioxide.l9 1.609e-03 1.939e-03 0.830 0.406663
## ammonia.l10 6.541e-02 1.369e-02 4.780 1.80e-06 ***
## hydrogensulfide.l10 -6.139e-03 1.275e-03 -4.814 1.51e-06 ***
## carbondioxide.l10 -8.184e-04 1.938e-03 -0.422 0.672916
## ammonia.l11 -1.302e-02 1.367e-02 -0.952 0.340973
## hydrogensulfide.l11 -1.044e-02 1.278e-03 -8.168 3.79e-16 ***
## carbondioxide.l11 3.511e-04 1.934e-03 0.182 0.855964
## ammonia.l12 8.622e-02 1.356e-02 6.360 2.17e-10 ***
## hydrogensulfide.l12 2.152e-03 1.281e-03 1.680 0.093093 .
## carbondioxide.l12 -1.214e-03 1.930e-03 -0.629 0.529430
## ammonia.l13 9.509e-03 1.350e-02 0.704 0.481257
## hydrogensulfide.l13 2.267e-03 1.280e-03 1.771 0.076545 .
## carbondioxide.l13 -2.774e-03 1.897e-03 -1.462 0.143740
## ammonia.l14 -3.028e-02 1.310e-02 -2.312 0.020838 *
## hydrogensulfide.l14 1.090e-02 1.275e-03 8.551 < 2e-16 ***
## carbondioxide.l14 -7.976e-04 1.608e-03 -0.496 0.619879
## const -4.713e-01 4.202e-01 -1.122 0.262094
## trend -3.624e-05 1.549e-05 -2.340 0.019316 *
## sd1 -1.520e-01 2.917e-01 -0.521 0.602303
## sd2 1.801e-01 2.922e-01 0.616 0.537620
## sd3 -1.639e-01 2.917e-01 -0.562 0.574242
## sd4 1.111e-01 2.918e-01 0.381 0.703319
## sd5 9.940e-02 2.918e-01 0.341 0.733405
## sd6 4.671e-01 2.918e-01 1.601 0.109485
## sd7 1.453e-01 2.915e-01 0.498 0.618234
## sd8 2.083e-01 2.920e-01 0.713 0.475601
## sd9 -4.776e-01 2.925e-01 -1.633 0.102566
## sd10 3.789e-01 2.920e-01 1.298 0.194464
## sd11 2.685e-01 2.915e-01 0.921 0.356969
## sd12 -1.073e-02 2.917e-01 -0.037 0.970663
## sd13 -3.162e-02 2.918e-01 -0.108 0.913722
## sd14 1.197e-01 2.915e-01 0.411 0.681420
## sd15 -2.079e-01 2.901e-01 -0.717 0.473620
## sd16 2.253e-01 2.904e-01 0.776 0.437902
## sd17 1.059e-01 2.905e-01 0.364 0.715516
## sd18 5.038e-02 2.899e-01 0.174 0.862048
## sd19 2.040e-01 2.910e-01 0.701 0.483313
## sd20 -1.007e-01 2.905e-01 -0.347 0.728875
## sd21 1.680e-01 2.904e-01 0.578 0.562968
## sd22 7.071e-01 2.901e-01 2.438 0.014807 *
## sd23 -1.139e-01 2.904e-01 -0.392 0.694851
## sd24 3.011e-01 2.902e-01 1.037 0.299575
## sd25 -6.850e-02 2.905e-01 -0.236 0.813593
## sd26 3.019e-01 2.899e-01 1.041 0.297705
## sd27 1.989e-01 2.904e-01 0.685 0.493343
## sd28 1.986e-02 2.900e-01 0.068 0.945422
## sd29 1.941e-01 2.898e-01 0.670 0.503154
## sd30 7.078e-01 2.904e-01 2.438 0.014810 *
## sd31 1.662e-01 2.901e-01 0.573 0.566676
## sd32 7.215e-02 2.914e-01 0.248 0.804481
## sd33 5.210e-02 2.917e-01 0.179 0.858267
## sd34 -1.170e-01 2.915e-01 -0.401 0.688078
## sd35 9.677e-02 2.913e-01 0.332 0.739734
## sd36 -1.847e-02 2.912e-01 -0.063 0.949441
## sd37 2.796e-01 2.911e-01 0.960 0.336866
## sd38 1.335e-01 2.913e-01 0.458 0.646860
## sd39 -2.801e-02 2.913e-01 -0.096 0.923393
## sd40 2.612e-01 2.908e-01 0.898 0.369272
## sd41 2.413e-01 2.915e-01 0.828 0.407824
## sd42 1.113e-01 2.911e-01 0.382 0.702213
## sd43 1.770e-02 2.911e-01 0.061 0.951534
## sd44 2.130e-01 2.914e-01 0.731 0.464990
## sd45 1.235e-01 2.912e-01 0.424 0.671385
## sd46 1.109e-01 2.910e-01 0.381 0.703222
## sd47 3.138e-01 2.910e-01 1.078 0.280971
## sd48 1.478e-01 2.913e-01 0.507 0.611859
## sd49 2.476e-01 2.910e-01 0.851 0.394780
## sd50 9.353e-02 2.915e-01 0.321 0.748337
## sd51 1.685e-01 2.911e-01 0.579 0.562808
## sd52 2.614e-01 2.913e-01 0.897 0.369679
## sd53 1.334e-01 2.910e-01 0.458 0.646707
## sd54 1.699e-01 2.912e-01 0.584 0.559548
## sd55 -3.782e-02 2.913e-01 -0.130 0.896708
## sd56 -4.665e-02 2.914e-01 -0.160 0.872804
## sd57 2.336e-01 2.910e-01 0.803 0.422039
## sd58 5.373e-01 2.913e-01 1.845 0.065141 .
## sd59 1.875e-01 2.910e-01 0.644 0.519370
## sd60 6.189e-01 2.911e-01 2.126 0.033563 *
## sd61 -2.298e-01 2.915e-01 -0.789 0.430384
## sd62 2.688e-01 2.913e-01 0.923 0.356172
## sd63 2.531e-01 2.914e-01 0.868 0.385222
## sd64 1.967e-01 2.913e-01 0.675 0.499587
## sd65 1.891e-01 2.915e-01 0.649 0.516520
## sd66 -2.345e-01 2.916e-01 -0.804 0.421313
## sd67 1.561e-01 2.913e-01 0.536 0.592107
## sd68 3.513e-02 2.913e-01 0.121 0.904039
## sd69 3.546e-02 2.911e-01 0.122 0.903054
## sd70 8.768e-02 2.913e-01 0.301 0.763478
## sd71 1.725e-01 2.915e-01 0.592 0.554017
## sd72 2.527e-01 2.913e-01 0.867 0.385715
## sd73 2.377e-01 2.914e-01 0.816 0.414649
## sd74 1.261e-01 2.914e-01 0.433 0.665055
## sd75 -9.159e-02 2.914e-01 -0.314 0.753306
## sd76 -2.732e-03 2.913e-01 -0.009 0.992517
## sd77 -3.340e-01 2.911e-01 -1.147 0.251245
## sd78 -9.867e-02 2.912e-01 -0.339 0.734741
## sd79 2.932e-01 2.909e-01 1.008 0.313580
## sd80 1.193e-01 2.911e-01 0.410 0.681798
## sd81 -1.812e-01 2.911e-01 -0.622 0.533705
## sd82 2.473e-01 2.913e-01 0.849 0.396051
## sd83 2.906e-01 2.911e-01 0.998 0.318288
## sd84 5.762e-01 2.911e-01 1.979 0.047817 *
## sd85 3.123e-01 2.912e-01 1.073 0.283501
## sd86 6.922e-01 2.916e-01 2.374 0.017648 *
## sd87 3.632e-01 2.916e-01 1.245 0.213102
## sd88 1.249e-01 2.913e-01 0.429 0.668107
## sd89 4.790e-01 2.912e-01 1.645 0.099997 .
## sd90 2.744e-01 2.915e-01 0.941 0.346644
## sd91 6.048e-01 2.912e-01 2.077 0.037859 *
## sd92 -1.608e-01 2.912e-01 -0.552 0.580957
## sd93 -3.942e-01 2.914e-01 -1.353 0.176189
## sd94 2.086e-01 2.918e-01 0.715 0.474641
## sd95 8.021e-02 2.917e-01 0.275 0.783333
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 1.626 on 5926 degrees of freedom
## Multiple R-Squared: 0.8439, Adjusted R-squared: 0.8402
## F-statistic: 232.1 on 138 and 5926 DF, p-value: < 2.2e-16
##
##
## Estimation results for equation hydrogensulfide:
## ================================================
## hydrogensulfide = ammonia.l1 + hydrogensulfide.l1 + carbondioxide.l1 + ammonia.l2 + hydrogensulfide.l2 + carbondioxide.l2 + ammonia.l3 + hydrogensulfide.l3 + carbondioxide.l3 + ammonia.l4 + hydrogensulfide.l4 + carbondioxide.l4 + ammonia.l5 + hydrogensulfide.l5 + carbondioxide.l5 + ammonia.l6 + hydrogensulfide.l6 + carbondioxide.l6 + ammonia.l7 + hydrogensulfide.l7 + carbondioxide.l7 + ammonia.l8 + hydrogensulfide.l8 + carbondioxide.l8 + ammonia.l9 + hydrogensulfide.l9 + carbondioxide.l9 + ammonia.l10 + hydrogensulfide.l10 + carbondioxide.l10 + ammonia.l11 + hydrogensulfide.l11 + carbondioxide.l11 + ammonia.l12 + hydrogensulfide.l12 + carbondioxide.l12 + ammonia.l13 + hydrogensulfide.l13 + carbondioxide.l13 + ammonia.l14 + hydrogensulfide.l14 + carbondioxide.l14 + const + trend + sd1 + sd2 + sd3 + sd4 + sd5 + sd6 + sd7 + sd8 + sd9 + sd10 + sd11 + sd12 + sd13 + sd14 + sd15 + sd16 + sd17 + sd18 + sd19 + sd20 + sd21 + sd22 + sd23 + sd24 + sd25 + sd26 + sd27 + sd28 + sd29 + sd30 + sd31 + sd32 + sd33 + sd34 + sd35 + sd36 + sd37 + sd38 + sd39 + sd40 + sd41 + sd42 + sd43 + sd44 + sd45 + sd46 + sd47 + sd48 + sd49 + sd50 + sd51 + sd52 + sd53 + sd54 + sd55 + sd56 + sd57 + sd58 + sd59 + sd60 + sd61 + sd62 + sd63 + sd64 + sd65 + sd66 + sd67 + sd68 + sd69 + sd70 + sd71 + sd72 + sd73 + sd74 + sd75 + sd76 + sd77 + sd78 + sd79 + sd80 + sd81 + sd82 + sd83 + sd84 + sd85 + sd86 + sd87 + sd88 + sd89 + sd90 + sd91 + sd92 + sd93 + sd94 + sd95
##
## Estimate Std. Error t value Pr(>|t|)
## ammonia.l1 1.659e-01 1.341e-01 1.237 0.216141
## hydrogensulfide.l1 9.544e-02 1.295e-02 7.372 1.91e-13 ***
## carbondioxide.l1 6.676e-03 1.643e-02 0.406 0.684515
## ammonia.l2 -7.970e-01 1.380e-01 -5.773 8.16e-09 ***
## hydrogensulfide.l2 -2.067e-02 1.300e-02 -1.589 0.112067
## carbondioxide.l2 -3.163e-02 1.948e-02 -1.624 0.104535
## ammonia.l3 2.184e-01 1.391e-01 1.571 0.116346
## hydrogensulfide.l3 1.115e-01 1.300e-02 8.578 < 2e-16 ***
## carbondioxide.l3 5.080e-02 1.980e-02 2.566 0.010323 *
## ammonia.l4 5.828e-02 1.396e-01 0.418 0.676237
## hydrogensulfide.l4 -5.197e-02 1.301e-02 -3.995 6.55e-05 ***
## carbondioxide.l4 3.789e-02 1.981e-02 1.913 0.055842 .
## ammonia.l5 4.380e-02 1.398e-01 0.313 0.754010
## hydrogensulfide.l5 3.826e-03 1.303e-02 0.294 0.769014
## carbondioxide.l5 -5.633e-02 1.979e-02 -2.846 0.004442 **
## ammonia.l6 7.430e-03 1.400e-01 0.053 0.957673
## hydrogensulfide.l6 -2.469e-02 1.303e-02 -1.895 0.058098 .
## carbondioxide.l6 7.799e-03 1.981e-02 0.394 0.693824
## ammonia.l7 3.515e-01 1.401e-01 2.509 0.012121 *
## hydrogensulfide.l7 1.802e-01 1.301e-02 13.845 < 2e-16 ***
## carbondioxide.l7 -2.070e-02 1.982e-02 -1.044 0.296365
## ammonia.l8 1.739e-01 1.401e-01 1.241 0.214495
## hydrogensulfide.l8 -5.459e-02 1.302e-02 -4.194 2.78e-05 ***
## carbondioxide.l8 2.091e-02 1.982e-02 1.055 0.291401
## ammonia.l9 5.229e-01 1.400e-01 3.736 0.000188 ***
## hydrogensulfide.l9 -1.516e-03 1.303e-02 -0.116 0.907348
## carbondioxide.l9 2.291e-02 1.980e-02 1.157 0.247297
## ammonia.l10 4.835e-01 1.398e-01 3.460 0.000545 ***
## hydrogensulfide.l10 -5.106e-02 1.302e-02 -3.921 8.92e-05 ***
## carbondioxide.l10 -3.681e-02 1.980e-02 -1.859 0.063058 .
## ammonia.l11 -4.675e-01 1.396e-01 -3.348 0.000819 ***
## hydrogensulfide.l11 9.509e-02 1.305e-02 7.288 3.56e-13 ***
## carbondioxide.l11 -1.049e-02 1.976e-02 -0.531 0.595599
## ammonia.l12 1.008e-01 1.385e-01 0.728 0.466772
## hydrogensulfide.l12 4.863e-02 1.308e-02 3.717 0.000204 ***
## carbondioxide.l12 1.899e-02 1.971e-02 0.963 0.335446
## ammonia.l13 -2.753e-01 1.379e-01 -1.997 0.045906 *
## hydrogensulfide.l13 -1.667e-02 1.307e-02 -1.275 0.202403
## carbondioxide.l13 5.683e-03 1.938e-02 0.293 0.769291
## ammonia.l14 -3.661e-01 1.338e-01 -2.736 0.006229 **
## hydrogensulfide.l14 1.374e-01 1.302e-02 10.549 < 2e-16 ***
## carbondioxide.l14 7.734e-03 1.642e-02 0.471 0.637674
## const -1.455e+01 4.292e+00 -3.390 0.000702 ***
## trend 1.807e-04 1.582e-04 1.143 0.253280
## sd1 -1.087e+00 2.979e+00 -0.365 0.715332
## sd2 -2.399e+00 2.984e+00 -0.804 0.421430
## sd3 -4.576e-01 2.979e+00 -0.154 0.877926
## sd4 -3.443e-01 2.980e+00 -0.116 0.908015
## sd5 2.893e+00 2.980e+00 0.971 0.331654
## sd6 -1.685e+00 2.980e+00 -0.566 0.571730
## sd7 -1.627e-01 2.978e+00 -0.055 0.956414
## sd8 8.148e-02 2.982e+00 0.027 0.978202
## sd9 2.679e+00 2.987e+00 0.897 0.369890
## sd10 -1.404e+00 2.982e+00 -0.471 0.637822
## sd11 -3.955e-01 2.977e+00 -0.133 0.894306
## sd12 -2.729e+00 2.979e+00 -0.916 0.359663
## sd13 7.810e+00 2.980e+00 2.620 0.008805 **
## sd14 -6.744e-01 2.977e+00 -0.227 0.820765
## sd15 9.623e-01 2.963e+00 0.325 0.745351
## sd16 8.074e-01 2.965e+00 0.272 0.785433
## sd17 1.585e+00 2.967e+00 0.534 0.593060
## sd18 2.051e+00 2.961e+00 0.693 0.488540
## sd19 4.975e-02 2.972e+00 0.017 0.986646
## sd20 -1.355e+00 2.967e+00 -0.457 0.647769
## sd21 3.484e-01 2.965e+00 0.117 0.906468
## sd22 1.751e-01 2.962e+00 0.059 0.952862
## sd23 9.140e-01 2.966e+00 0.308 0.757955
## sd24 1.026e+00 2.964e+00 0.346 0.729228
## sd25 -3.838e-01 2.967e+00 -0.129 0.897068
## sd26 6.353e-01 2.961e+00 0.215 0.830124
## sd27 1.489e+00 2.965e+00 0.502 0.615699
## sd28 -5.252e-01 2.962e+00 -0.177 0.859272
## sd29 5.906e-01 2.960e+00 0.200 0.841868
## sd30 -1.523e-01 2.965e+00 -0.051 0.959050
## sd31 -3.943e-01 2.963e+00 -0.133 0.894128
## sd32 5.506e-01 2.976e+00 0.185 0.853247
## sd33 1.093e+00 2.980e+00 0.367 0.713868
## sd34 -9.401e-01 2.977e+00 -0.316 0.752169
## sd35 9.802e-01 2.975e+00 0.329 0.741806
## sd36 -3.678e-01 2.974e+00 -0.124 0.901588
## sd37 -4.965e-01 2.973e+00 -0.167 0.867391
## sd38 4.515e+00 2.975e+00 1.518 0.129159
## sd39 -7.837e-01 2.975e+00 -0.263 0.792193
## sd40 1.101e+00 2.970e+00 0.371 0.711011
## sd41 -1.177e+00 2.977e+00 -0.395 0.692691
## sd42 1.733e+00 2.973e+00 0.583 0.559996
## sd43 -1.477e-02 2.973e+00 -0.005 0.996036
## sd44 5.841e+00 2.977e+00 1.962 0.049786 *
## sd45 -1.078e+00 2.974e+00 -0.362 0.717003
## sd46 4.384e-01 2.972e+00 0.147 0.882765
## sd47 -1.231e+00 2.972e+00 -0.414 0.678696
## sd48 4.675e-01 2.975e+00 0.157 0.875156
## sd49 -7.899e-01 2.972e+00 -0.266 0.790402
## sd50 7.144e-01 2.977e+00 0.240 0.810360
## sd51 -9.325e-01 2.973e+00 -0.314 0.753813
## sd52 -2.741e-01 2.975e+00 -0.092 0.926599
## sd53 -1.663e-01 2.972e+00 -0.056 0.955388
## sd54 2.421e-01 2.974e+00 0.081 0.935106
## sd55 3.155e-01 2.975e+00 0.106 0.915531
## sd56 -3.814e-01 2.976e+00 -0.128 0.898026
## sd57 -1.661e-01 2.972e+00 -0.056 0.955425
## sd58 2.501e+00 2.975e+00 0.841 0.400492
## sd59 -4.874e-01 2.972e+00 -0.164 0.869757
## sd60 1.360e+00 2.973e+00 0.457 0.647448
## sd61 -5.622e-01 2.977e+00 -0.189 0.850202
## sd62 6.130e-01 2.975e+00 0.206 0.836778
## sd63 -3.109e-01 2.976e+00 -0.104 0.916802
## sd64 6.895e-01 2.976e+00 0.232 0.816750
## sd65 4.459e-02 2.977e+00 0.015 0.988049
## sd66 6.555e-01 2.979e+00 0.220 0.825815
## sd67 -6.332e-01 2.975e+00 -0.213 0.831434
## sd68 1.010e+00 2.976e+00 0.339 0.734322
## sd69 -7.061e-01 2.973e+00 -0.238 0.812261
## sd70 -1.184e-01 2.976e+00 -0.040 0.968252
## sd71 -3.927e-01 2.977e+00 -0.132 0.895044
## sd72 -4.032e-01 2.975e+00 -0.136 0.892201
## sd73 -6.422e-02 2.976e+00 -0.022 0.982784
## sd74 8.019e-01 2.976e+00 0.269 0.787567
## sd75 9.609e-02 2.976e+00 0.032 0.974244
## sd76 3.824e-01 2.975e+00 0.129 0.897743
## sd77 -1.769e-01 2.973e+00 -0.059 0.952559
## sd78 1.665e+00 2.974e+00 0.560 0.575521
## sd79 -7.014e-01 2.971e+00 -0.236 0.813406
## sd80 -6.638e-01 2.973e+00 -0.223 0.823294
## sd81 -2.071e-01 2.973e+00 -0.070 0.944455
## sd82 -2.312e-01 2.975e+00 -0.078 0.938072
## sd83 5.051e+00 2.973e+00 1.699 0.089430 .
## sd84 1.131e-01 2.973e+00 0.038 0.969652
## sd85 1.203e-01 2.974e+00 0.040 0.967737
## sd86 2.131e+00 2.979e+00 0.715 0.474363
## sd87 9.235e+00 2.979e+00 3.100 0.001942 **
## sd88 4.496e+00 2.975e+00 1.511 0.130823
## sd89 1.261e+00 2.974e+00 0.424 0.671625
## sd90 -2.065e+00 2.977e+00 -0.694 0.487913
## sd91 9.099e+00 2.974e+00 3.060 0.002227 **
## sd92 1.048e-01 2.974e+00 0.035 0.971884
## sd93 1.110e+00 2.976e+00 0.373 0.709180
## sd94 1.136e+01 2.980e+00 3.812 0.000139 ***
## sd95 4.994e+00 2.979e+00 1.677 0.093686 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 16.6 on 5926 degrees of freedom
## Multiple R-Squared: 0.1379, Adjusted R-squared: 0.1178
## F-statistic: 6.867 on 138 and 5926 DF, p-value: < 2.2e-16
##
##
## Estimation results for equation carbondioxide:
## ==============================================
## carbondioxide = ammonia.l1 + hydrogensulfide.l1 + carbondioxide.l1 + ammonia.l2 + hydrogensulfide.l2 + carbondioxide.l2 + ammonia.l3 + hydrogensulfide.l3 + carbondioxide.l3 + ammonia.l4 + hydrogensulfide.l4 + carbondioxide.l4 + ammonia.l5 + hydrogensulfide.l5 + carbondioxide.l5 + ammonia.l6 + hydrogensulfide.l6 + carbondioxide.l6 + ammonia.l7 + hydrogensulfide.l7 + carbondioxide.l7 + ammonia.l8 + hydrogensulfide.l8 + carbondioxide.l8 + ammonia.l9 + hydrogensulfide.l9 + carbondioxide.l9 + ammonia.l10 + hydrogensulfide.l10 + carbondioxide.l10 + ammonia.l11 + hydrogensulfide.l11 + carbondioxide.l11 + ammonia.l12 + hydrogensulfide.l12 + carbondioxide.l12 + ammonia.l13 + hydrogensulfide.l13 + carbondioxide.l13 + ammonia.l14 + hydrogensulfide.l14 + carbondioxide.l14 + const + trend + sd1 + sd2 + sd3 + sd4 + sd5 + sd6 + sd7 + sd8 + sd9 + sd10 + sd11 + sd12 + sd13 + sd14 + sd15 + sd16 + sd17 + sd18 + sd19 + sd20 + sd21 + sd22 + sd23 + sd24 + sd25 + sd26 + sd27 + sd28 + sd29 + sd30 + sd31 + sd32 + sd33 + sd34 + sd35 + sd36 + sd37 + sd38 + sd39 + sd40 + sd41 + sd42 + sd43 + sd44 + sd45 + sd46 + sd47 + sd48 + sd49 + sd50 + sd51 + sd52 + sd53 + sd54 + sd55 + sd56 + sd57 + sd58 + sd59 + sd60 + sd61 + sd62 + sd63 + sd64 + sd65 + sd66 + sd67 + sd68 + sd69 + sd70 + sd71 + sd72 + sd73 + sd74 + sd75 + sd76 + sd77 + sd78 + sd79 + sd80 + sd81 + sd82 + sd83 + sd84 + sd85 + sd86 + sd87 + sd88 + sd89 + sd90 + sd91 + sd92 + sd93 + sd94 + sd95
##
## Estimate Std. Error t value Pr(>|t|)
## ammonia.l1 -0.4234275 0.1072126 -3.949 7.92e-05 ***
## hydrogensulfide.l1 0.0180098 0.0103472 1.741 0.081815 .
## carbondioxide.l1 0.6438456 0.0131321 49.029 < 2e-16 ***
## ammonia.l2 0.1177476 0.1103258 1.067 0.285892
## hydrogensulfide.l2 -0.0194226 0.0103937 -1.869 0.061714 .
## carbondioxide.l2 0.2194782 0.0155711 14.095 < 2e-16 ***
## ammonia.l3 -0.2266287 0.1111597 -2.039 0.041517 *
## hydrogensulfide.l3 -0.0267058 0.0103934 -2.570 0.010209 *
## carbondioxide.l3 0.0184779 0.0158243 1.168 0.242977
## ammonia.l4 0.4336624 0.1115426 3.888 0.000102 ***
## hydrogensulfide.l4 -0.0396459 0.0103975 -3.813 0.000139 ***
## carbondioxide.l4 -0.0062392 0.0158319 -0.394 0.693528
## ammonia.l5 0.0644197 0.1117199 0.577 0.564219
## hydrogensulfide.l5 0.0253858 0.0104128 2.438 0.014801 *
## carbondioxide.l5 0.0330747 0.0158198 2.091 0.036595 *
## ammonia.l6 0.2105700 0.1118879 1.882 0.059889 .
## hydrogensulfide.l6 0.0197087 0.0104135 1.893 0.058458 .
## carbondioxide.l6 0.0505471 0.0158333 3.192 0.001418 **
## ammonia.l7 0.0679006 0.1119367 0.607 0.544141
## hydrogensulfide.l7 0.0055220 0.0104004 0.531 0.595479
## carbondioxide.l7 -0.0114000 0.0158400 -0.720 0.471742
## ammonia.l8 -0.1941234 0.1119715 -1.734 0.083026 .
## hydrogensulfide.l8 0.0004894 0.0104021 0.047 0.962473
## carbondioxide.l8 0.0033246 0.0158374 0.210 0.833737
## ammonia.l9 0.1165815 0.1118545 1.042 0.297333
## hydrogensulfide.l9 -0.0239031 0.0104100 -2.296 0.021701 *
## carbondioxide.l9 -0.0316721 0.0158265 -2.001 0.045415 *
## ammonia.l10 -0.1272829 0.1117085 -1.139 0.254574
## hydrogensulfide.l10 0.0159606 0.0104075 1.534 0.125189
## carbondioxide.l10 0.0167994 0.0158226 1.062 0.288400
## ammonia.l11 0.1317935 0.1115938 1.181 0.237646
## hydrogensulfide.l11 -0.0353883 0.0104285 -3.393 0.000695 ***
## carbondioxide.l11 -0.0007360 0.0157904 -0.047 0.962826
## ammonia.l12 -0.1550349 0.1106556 -1.401 0.161249
## hydrogensulfide.l12 0.0088599 0.0104568 0.847 0.396870
## carbondioxide.l12 -0.0197873 0.0157526 -1.256 0.209118
## ammonia.l13 -0.0031826 0.1102072 -0.029 0.976963
## hydrogensulfide.l13 0.0021914 0.0104478 0.210 0.833869
## carbondioxide.l13 0.0018226 0.0154862 0.118 0.906316
## ammonia.l14 0.0769552 0.1069353 0.720 0.471773
## hydrogensulfide.l14 -0.0021766 0.0104063 -0.209 0.834328
## carbondioxide.l14 0.0326431 0.0131249 2.487 0.012906 *
## const 28.1148962 3.4302246 8.196 3.01e-16 ***
## trend -0.0002768 0.0001264 -2.189 0.028622 *
## sd1 -0.5379689 2.3810925 -0.226 0.821261
## sd2 -1.2443674 2.3848016 -0.522 0.601836
## sd3 0.0818930 2.3812042 0.034 0.972566
## sd4 0.7951990 2.3816243 0.334 0.738475
## sd5 1.9159915 2.3819250 0.804 0.421205
## sd6 0.6128484 2.3819140 0.257 0.796962
## sd7 -1.8540013 2.3797897 -0.779 0.435975
## sd8 -4.7419293 2.3831491 -1.990 0.046662 *
## sd9 -2.8980343 2.3876045 -1.214 0.224879
## sd10 0.4738645 2.3831426 0.199 0.842395
## sd11 0.9442525 2.3792556 0.397 0.691478
## sd12 -2.6766374 2.3806906 -1.124 0.260927
## sd13 0.3369924 2.3820706 0.141 0.887503
## sd14 1.4074953 2.3791130 0.592 0.554138
## sd15 0.1616304 2.3679496 0.068 0.945583
## sd16 0.9767906 2.3700997 0.412 0.680259
## sd17 -0.4927691 2.3710376 -0.208 0.835370
## sd18 0.0424737 2.3663867 0.018 0.985680
## sd19 1.6820537 2.3753136 0.708 0.478886
## sd20 -1.1415073 2.3710109 -0.481 0.630219
## sd21 -2.7794740 2.3700407 -1.173 0.240942
## sd22 -2.5330443 2.3676857 -1.070 0.284735
## sd23 -2.6103502 2.3704682 -1.101 0.270856
## sd24 1.9879780 2.3691748 0.839 0.401446
## sd25 -0.2985414 2.3712864 -0.126 0.899817
## sd26 -0.0889686 2.3664099 -0.038 0.970011
## sd27 -0.3756527 2.3700639 -0.158 0.874069
## sd28 -0.1840984 2.3673773 -0.078 0.938018
## sd29 0.7733254 2.3658317 0.327 0.743776
## sd30 0.0222857 2.3700808 0.009 0.992498
## sd31 3.9391696 2.3678732 1.664 0.096247 .
## sd32 -0.7783652 2.3789001 -0.327 0.743532
## sd33 0.2624532 2.3813804 0.110 0.912246
## sd34 0.4495148 2.3791423 0.189 0.850146
## sd35 -3.1476264 2.3776212 -1.324 0.185602
## sd36 -3.4336472 2.3769350 -1.445 0.148632
## sd37 0.5795289 2.3764469 0.244 0.807345
## sd38 0.1035992 2.3777675 0.044 0.965249
## sd39 -2.4017488 2.3774319 -1.010 0.312427
## sd40 -1.1823128 2.3740611 -0.498 0.618494
## sd41 -1.6570666 2.3794662 -0.696 0.486204
## sd42 -0.1086294 2.3757714 -0.046 0.963532
## sd43 0.1719526 2.3764040 0.072 0.942319
## sd44 1.3749047 2.3789626 0.578 0.563325
## sd45 -0.3638319 2.3768303 -0.153 0.878345
## sd46 -1.3830584 2.3756875 -0.582 0.560473
## sd47 -1.3979353 2.3756648 -0.588 0.556260
## sd48 0.6347112 2.3779510 0.267 0.789544
## sd49 1.0921006 2.3752584 0.460 0.645690
## sd50 -0.8833846 2.3793473 -0.371 0.710448
## sd51 1.6333178 2.3762490 0.687 0.491888
## sd52 0.4360420 2.3780835 0.183 0.854523
## sd53 -3.2064487 2.3751911 -1.350 0.177076
## sd54 -1.1965825 2.3767368 -0.503 0.614662
## sd55 2.0460748 2.3776009 0.861 0.389514
## sd56 -2.2113002 2.3783987 -0.930 0.352542
## sd57 0.5331097 2.3751134 0.224 0.822410
## sd58 0.0975357 2.3775107 0.041 0.967278
## sd59 -1.6654063 2.3756040 -0.701 0.483302
## sd60 0.1649883 2.3764832 0.069 0.944653
## sd61 -0.4548878 2.3790027 -0.191 0.848368
## sd62 5.1004809 2.3779144 2.145 0.031998 *
## sd63 1.1175237 2.3788970 0.470 0.638540
## sd64 -1.6648732 2.3781628 -0.700 0.483913
## sd65 -0.4418780 2.3791077 -0.186 0.852661
## sd66 -1.5925550 2.3805220 -0.669 0.503525
## sd67 0.3804215 2.3774330 0.160 0.872876
## sd68 0.6508273 2.3781443 0.274 0.784348
## sd69 -0.5606870 2.3760753 -0.236 0.813463
## sd70 2.6298598 2.3781627 1.106 0.268842
## sd71 0.0557424 2.3792318 0.023 0.981309
## sd72 1.2388323 2.3777838 0.521 0.602384
## sd73 -0.3826507 2.3785286 -0.161 0.872196
## sd74 -0.7253846 2.3781669 -0.305 0.760363
## sd75 -2.9791587 2.3786087 -1.252 0.210445
## sd76 -0.4860843 2.3779870 -0.204 0.838040
## sd77 -0.2647202 2.3761603 -0.111 0.911298
## sd78 -0.0201609 2.3768197 -0.008 0.993232
## sd79 -0.0564806 2.3747788 -0.024 0.981026
## sd80 -1.0167303 2.3757905 -0.428 0.668700
## sd81 -1.0627282 2.3759203 -0.447 0.654681
## sd82 0.6478705 2.3781100 0.272 0.785300
## sd83 -0.3547049 2.3764847 -0.149 0.881357
## sd84 -0.1000723 2.3762749 -0.042 0.966410
## sd85 2.8107567 2.3767048 1.183 0.237004
## sd86 0.7779452 2.3805505 0.327 0.743837
## sd87 1.1315008 2.3805778 0.475 0.634587
## sd88 -3.4025042 2.3778876 -1.431 0.152513
## sd89 -0.2836936 2.3767311 -0.119 0.904992
## sd90 0.8073528 2.3795283 0.339 0.734402
## sd91 -0.3690178 2.3769890 -0.155 0.876633
## sd92 -1.8389141 2.3771822 -0.774 0.439217
## sd93 -3.9284278 2.3785367 -1.652 0.098666 .
## sd94 0.5565526 2.3818475 0.234 0.815254
## sd95 -0.7672144 2.3807786 -0.322 0.747272
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 13.27 on 5926 degrees of freedom
## Multiple R-Squared: 0.8722, Adjusted R-squared: 0.8692
## F-statistic: 293.1 on 138 and 5926 DF, p-value: < 2.2e-16
##
##
##
## Covariance matrix of residuals:
## ammonia hydrogensulfide carbondioxide
## ammonia 2.642 2.5613 3.2642
## hydrogensulfide 2.561 275.6101 -0.7422
## carbondioxide 3.264 -0.7422 176.0517
##
## Correlation matrix of residuals:
## ammonia hydrogensulfide carbondioxide
## ammonia 1.00000 0.09491 0.15134
## hydrogensulfide 0.09491 1.00000 -0.00337
## carbondioxide 0.15134 -0.00337 1.00000var_res_fit <- restrict(var_fit, method = "ser", thresh = 2)
summary(var_res_fit)##
## VAR Estimation Results:
## =========================
## Endogenous variables: ammonia, hydrogensulfide, carbondioxide
## Deterministic variables: both
## Sample size: 6065
## Log Likelihood: -61280.729
## Roots of the characteristic polynomial:
## 0.9912 0.9714 0.9254 0.9013 0.9013 0.8998 0.8998 0.8977 0.8977 0.8591 0.8591 0.8495 0.8495 0.8298 0.8298 0.8211 0.8211 0.8111 0.8111 0.8064 0.8064 0.799 0.799 0.798 0.798 0.7939 0.7939 0.7687 0.7687 0.7569 0.7559 0.7559 0.7546 0.7546 0.7536 0.7536 0.7412 0.7412 0.7355 0.7355 0.1556 0.06387
## Call:
## VAR(y = gas, p = 14, type = "both", season = 96L)
##
##
## Estimation results for equation ammonia:
## ========================================
## ammonia = ammonia.l1 + hydrogensulfide.l1 + carbondioxide.l1 + ammonia.l2 + hydrogensulfide.l2 + ammonia.l3 + carbondioxide.l3 + ammonia.l4 + ammonia.l5 + carbondioxide.l5 + ammonia.l6 + ammonia.l8 + carbondioxide.l8 + ammonia.l9 + ammonia.l10 + hydrogensulfide.l10 + hydrogensulfide.l11 + ammonia.l12 + hydrogensulfide.l13 + carbondioxide.l13 + ammonia.l14 + hydrogensulfide.l14 + trend + sd9 + sd22 + sd30 + sd60 + sd77 + sd84 + sd86 + sd91 + sd93
##
## Estimate Std. Error t value Pr(>|t|)
## ammonia.l1 2.296e-01 1.276e-02 17.988 < 2e-16 ***
## hydrogensulfide.l1 4.336e-03 1.206e-03 3.595 0.000327 ***
## carbondioxide.l1 1.006e-02 1.307e-03 7.693 1.66e-14 ***
## ammonia.l2 1.299e-01 1.313e-02 9.894 < 2e-16 ***
## hydrogensulfide.l2 -5.339e-03 1.212e-03 -4.404 1.08e-05 ***
## ammonia.l3 1.519e-01 1.322e-02 11.491 < 2e-16 ***
## carbondioxide.l3 -3.430e-03 1.618e-03 -2.120 0.034059 *
## ammonia.l4 1.002e-01 1.316e-02 7.613 3.09e-14 ***
## ammonia.l5 8.132e-02 1.325e-02 6.135 9.03e-10 ***
## carbondioxide.l5 -4.671e-03 1.520e-03 -3.072 0.002136 **
## ammonia.l6 6.249e-02 1.298e-02 4.815 1.51e-06 ***
## ammonia.l8 3.281e-02 1.312e-02 2.500 0.012458 *
## carbondioxide.l8 3.019e-03 1.337e-03 2.257 0.024021 *
## ammonia.l9 5.225e-02 1.311e-02 3.985 6.82e-05 ***
## ammonia.l10 6.429e-02 1.286e-02 4.998 5.96e-07 ***
## hydrogensulfide.l10 -5.862e-03 1.213e-03 -4.834 1.37e-06 ***
## hydrogensulfide.l11 -1.034e-02 1.208e-03 -8.559 < 2e-16 ***
## ammonia.l12 8.828e-02 1.253e-02 7.047 2.04e-12 ***
## hydrogensulfide.l13 2.902e-03 1.219e-03 2.380 0.017341 *
## carbondioxide.l13 -4.257e-03 9.946e-04 -4.280 1.90e-05 ***
## ammonia.l14 -2.970e-02 1.212e-02 -2.451 0.014285 *
## hydrogensulfide.l14 1.115e-02 1.219e-03 9.152 < 2e-16 ***
## trend -4.246e-05 1.459e-05 -2.911 0.003620 **
## sd9 -6.092e-01 2.068e-01 -2.946 0.003227 **
## sd22 5.671e-01 2.040e-01 2.780 0.005456 **
## sd30 5.777e-01 2.040e-01 2.831 0.004650 **
## sd60 4.984e-01 2.056e-01 2.424 0.015373 *
## sd77 -4.645e-01 2.055e-01 -2.260 0.023827 *
## sd84 4.350e-01 2.056e-01 2.115 0.034433 *
## sd86 5.733e-01 2.057e-01 2.786 0.005345 **
## sd91 4.821e-01 2.058e-01 2.343 0.019185 *
## sd93 -5.091e-01 2.060e-01 -2.471 0.013500 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 1.621 on 6033 degrees of freedom
## Multiple R-Squared: 0.9632, Adjusted R-squared: 0.963
## F-statistic: 4928 on 32 and 6033 DF, p-value: < 2.2e-16
##
##
## Estimation results for equation hydrogensulfide:
## ================================================
## hydrogensulfide = hydrogensulfide.l1 + ammonia.l2 + ammonia.l3 + hydrogensulfide.l3 + carbondioxide.l3 + hydrogensulfide.l4 + carbondioxide.l4 + carbondioxide.l5 + hydrogensulfide.l6 + ammonia.l7 + hydrogensulfide.l7 + hydrogensulfide.l8 + ammonia.l9 + ammonia.l10 + hydrogensulfide.l10 + ammonia.l11 + hydrogensulfide.l11 + hydrogensulfide.l12 + ammonia.l14 + hydrogensulfide.l14 + const + sd13 + sd38 + sd44 + sd83 + sd87 + sd88 + sd91 + sd94 + sd95
##
## Estimate Std. Error t value Pr(>|t|)
## hydrogensulfide.l1 0.09332 0.01257 7.422 1.31e-13 ***
## ammonia.l2 -0.75567 0.12452 -6.069 1.37e-09 ***
## ammonia.l3 0.30147 0.12716 2.371 0.017778 *
## hydrogensulfide.l3 0.10507 0.01261 8.331 < 2e-16 ***
## carbondioxide.l3 0.03634 0.01599 2.272 0.023127 *
## hydrogensulfide.l4 -0.05034 0.01263 -3.987 6.77e-05 ***
## carbondioxide.l4 0.03741 0.01853 2.018 0.043601 *
## carbondioxide.l5 -0.05314 0.01591 -3.341 0.000841 ***
## hydrogensulfide.l6 -0.02647 0.01224 -2.162 0.030643 *
## ammonia.l7 0.38988 0.12678 3.075 0.002113 **
## hydrogensulfide.l7 0.17813 0.01272 14.001 < 2e-16 ***
## hydrogensulfide.l8 -0.05121 0.01254 -4.084 4.49e-05 ***
## ammonia.l9 0.54946 0.12907 4.257 2.10e-05 ***
## ammonia.l10 0.47697 0.13201 3.613 0.000305 ***
## hydrogensulfide.l10 -0.05000 0.01252 -3.995 6.55e-05 ***
## ammonia.l11 -0.44126 0.12898 -3.421 0.000628 ***
## hydrogensulfide.l11 0.09669 0.01269 7.621 2.90e-14 ***
## hydrogensulfide.l12 0.04503 0.01227 3.671 0.000244 ***
## ammonia.l14 -0.35616 0.12100 -2.943 0.003259 **
## hydrogensulfide.l14 0.13746 0.01265 10.870 < 2e-16 ***
## const -12.03737 3.67584 -3.275 0.001064 **
## sd13 7.70652 2.09747 3.674 0.000241 ***
## sd38 4.64414 2.09379 2.218 0.026588 *
## sd44 5.56117 2.09474 2.655 0.007956 **
## sd83 4.96650 2.09427 2.371 0.017748 *
## sd87 9.17705 2.09584 4.379 1.21e-05 ***
## sd88 4.51068 2.09910 2.149 0.031684 *
## sd91 9.21995 2.10034 4.390 1.15e-05 ***
## sd94 11.54317 2.10041 5.496 4.05e-08 ***
## sd95 5.19383 2.10579 2.466 0.013673 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 16.51 on 6035 degrees of freedom
## Multiple R-Squared: 0.1364, Adjusted R-squared: 0.1321
## F-statistic: 31.76 on 30 and 6035 DF, p-value: < 2.2e-16
##
##
## Estimation results for equation carbondioxide:
## ==============================================
## carbondioxide = ammonia.l1 + carbondioxide.l1 + carbondioxide.l2 + hydrogensulfide.l3 + ammonia.l4 + hydrogensulfide.l4 + hydrogensulfide.l5 + carbondioxide.l5 + hydrogensulfide.l6 + carbondioxide.l6 + hydrogensulfide.l9 + carbondioxide.l9 + hydrogensulfide.l11 + carbondioxide.l14 + const + trend + sd8 + sd31 + sd62 + sd93
##
## Estimate Std. Error t value Pr(>|t|)
## ammonia.l1 -0.3747870 0.0842795 -4.447 8.87e-06 ***
## carbondioxide.l1 0.6433491 0.0126400 50.898 < 2e-16 ***
## carbondioxide.l2 0.2262035 0.0135347 16.713 < 2e-16 ***
## hydrogensulfide.l3 -0.0294291 0.0097989 -3.003 0.002682 **
## ammonia.l4 0.4644258 0.0845331 5.494 4.09e-08 ***
## hydrogensulfide.l4 -0.0412353 0.0099774 -4.133 3.63e-05 ***
## hydrogensulfide.l5 0.0252515 0.0098054 2.575 0.010040 *
## carbondioxide.l5 0.0352969 0.0136319 2.589 0.009640 **
## hydrogensulfide.l6 0.0256280 0.0098888 2.592 0.009576 **
## carbondioxide.l6 0.0454284 0.0136079 3.338 0.000848 ***
## hydrogensulfide.l9 -0.0283851 0.0097786 -2.903 0.003712 **
## carbondioxide.l9 -0.0273858 0.0107445 -2.549 0.010833 *
## hydrogensulfide.l11 -0.0315452 0.0098757 -3.194 0.001409 **
## carbondioxide.l14 0.0284592 0.0082714 3.441 0.000584 ***
## const 27.5144940 3.3889415 8.119 5.65e-16 ***
## trend -0.0002819 0.0001232 -2.289 0.022139 *
## sd8 -4.1713571 1.6798446 -2.483 0.013048 *
## sd31 4.1848904 1.6649206 2.514 0.011977 *
## sd62 5.3289989 1.6782201 3.175 0.001504 **
## sd93 -3.5281899 1.6800403 -2.100 0.035765 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 13.24 on 6045 degrees of freedom
## Multiple R-Squared: 0.9994, Adjusted R-squared: 0.9994
## F-statistic: 5.445e+05 on 20 and 6045 DF, p-value: < 2.2e-16
##
##
##
## Covariance matrix of residuals:
## ammonia hydrogensulfide carbondioxide
## ammonia 2.674 2.5519 3.3343
## hydrogensulfide 2.552 277.7304 -0.7236
## carbondioxide 3.334 -0.7236 178.7849
##
## Correlation matrix of residuals:
## ammonia hydrogensulfide carbondioxide
## ammonia 1.00000 0.093634 0.152484
## hydrogensulfide 0.09363 1.000000 -0.003248
## carbondioxide 0.15248 -0.003248 1.000000fit_ts <- ts(fitted(var_res_fit)[,1], frequency = frequency(ammonia_ts))
rmse <- accuracy(ammonia_ts, fit_ts)[,"RMSE"]
plot(ammonia_ts, ylab = 'ammonia', lwd = 1, lty ='dashed'); grid(); lines(fit_ts, col = 'red', lty = 'dashed')
legend("topright", c("ammonia", "fitted"), fill=c("black", "red")); mape <- accuracy(ammonia_ts, fit_ts)[,"RMSE"]
mes <- paste("RMSE = ", round(mape, 2), sep=''); text(60, 30, mes)
fit_ts <- ts(fitted(var_res_fit)[,2], frequency = frequency(hydrogensulfide_ts))
rmse <- accuracy(hydrogensulfide_ts, fit_ts)[,"RMSE"]
plot(hydrogensulfide_ts, ylab = 'hydrogensulfide', lwd = 1, lty ='dashed'); grid(); lines(fit_ts, col = 'red', lty = 'dashed')
legend("topright", c("hydrogensulfide", "fitted"), fill=c("black", "red")); mape <- accuracy(hydrogensulfide_ts, fit_ts)[,"RMSE"]
mes <- paste("RMSE = ", round(mape, 2), sep=''); text(60, 200, mes)
fit_ts <- ts(fitted(var_res_fit)[,3], frequency = frequency(carbondioxide_ts))
rmse <- accuracy(carbondioxide_ts, fit_ts)[,"RMSE"]
plot(carbondioxide_ts, ylab = 'carbondioxide', lwd = 1, lty ='dashed'); grid(); lines(fit_ts, col = 'red', lty = 'dashed')
legend("topright", c("carbondioxide", "fitted"), fill=c("black", "red")); mape <- accuracy(carbondioxide_ts, fit_ts)[,"RMSE"]
mes <- paste("RMSE = ", round(mape, 2), sep=''); text(60, 700, mes)
Acoef(var_res_fit)## [[1]]
## ammonia.l1 hydrogensulfide.l1 carbondioxide.l1
## ammonia 0.2296063 0.004336139 0.01005772
## hydrogensulfide 0.0000000 0.093318891 0.00000000
## carbondioxide -0.3747870 0.000000000 0.64334906
##
## [[2]]
## ammonia.l2 hydrogensulfide.l2 carbondioxide.l2
## ammonia 0.1299163 -0.005338731 0.0000000
## hydrogensulfide -0.7556705 0.000000000 0.0000000
## carbondioxide 0.0000000 0.000000000 0.2262035
##
## [[3]]
## ammonia.l3 hydrogensulfide.l3 carbondioxide.l3
## ammonia 0.1518805 0.00000000 -0.003429757
## hydrogensulfide 0.3014658 0.10507247 0.036335905
## carbondioxide 0.0000000 -0.02942909 0.000000000
##
## [[4]]
## ammonia.l4 hydrogensulfide.l4 carbondioxide.l4
## ammonia 0.1002091 0.00000000 0.00000000
## hydrogensulfide 0.0000000 -0.05033974 0.03740788
## carbondioxide 0.4644258 -0.04123528 0.00000000
##
## [[5]]
## ammonia.l5 hydrogensulfide.l5 carbondioxide.l5
## ammonia 0.08131656 0.00000000 -0.004670779
## hydrogensulfide 0.00000000 0.00000000 -0.053142094
## carbondioxide 0.00000000 0.02525148 0.035296948
##
## [[6]]
## ammonia.l6 hydrogensulfide.l6 carbondioxide.l6
## ammonia 0.06248569 0.00000000 0.00000000
## hydrogensulfide 0.00000000 -0.02646683 0.00000000
## carbondioxide 0.00000000 0.02562801 0.04542839
##
## [[7]]
## ammonia.l7 hydrogensulfide.l7 carbondioxide.l7
## ammonia 0.0000000 0.0000000 0
## hydrogensulfide 0.3898805 0.1781276 0
## carbondioxide 0.0000000 0.0000000 0
##
## [[8]]
## ammonia.l8 hydrogensulfide.l8 carbondioxide.l8
## ammonia 0.03280677 0.0000000 0.003018706
## hydrogensulfide 0.00000000 -0.0512121 0.000000000
## carbondioxide 0.00000000 0.0000000 0.000000000
##
## [[9]]
## ammonia.l9 hydrogensulfide.l9 carbondioxide.l9
## ammonia 0.0522469 0.0000000 0.0000000
## hydrogensulfide 0.5494619 0.0000000 0.0000000
## carbondioxide 0.0000000 -0.0283851 -0.0273858
##
## [[10]]
## ammonia.l10 hydrogensulfide.l10 carbondioxide.l10
## ammonia 0.06429442 -0.005861954 0
## hydrogensulfide 0.47696749 -0.050001370 0
## carbondioxide 0.00000000 0.000000000 0
##
## [[11]]
## ammonia.l11 hydrogensulfide.l11 carbondioxide.l11
## ammonia 0.0000000 -0.01033689 0
## hydrogensulfide -0.4412632 0.09668784 0
## carbondioxide 0.0000000 -0.03154515 0
##
## [[12]]
## ammonia.l12 hydrogensulfide.l12 carbondioxide.l12
## ammonia 0.08828156 0.0000000 0
## hydrogensulfide 0.00000000 0.0450265 0
## carbondioxide 0.00000000 0.0000000 0
##
## [[13]]
## ammonia.l13 hydrogensulfide.l13 carbondioxide.l13
## ammonia 0 0.002902069 -0.004256983
## hydrogensulfide 0 0.000000000 0.000000000
## carbondioxide 0 0.000000000 0.000000000
##
## [[14]]
## ammonia.l14 hydrogensulfide.l14 carbondioxide.l14
## ammonia -0.02969621 0.01115296 0.00000000
## hydrogensulfide -0.35615903 0.13745783 0.00000000
## carbondioxide 0.00000000 0.00000000 0.02845924roots(var_res_fit)## [1] 0.99123496 0.97139797 0.92540183 0.90132113 0.90132113 0.89978067
## [7] 0.89978067 0.89771295 0.89771295 0.85906015 0.85906015 0.84947036
## [13] 0.84947036 0.82975641 0.82975641 0.82112993 0.82112993 0.81107116
## [19] 0.81107116 0.80643102 0.80643102 0.79895743 0.79895743 0.79795516
## [25] 0.79795516 0.79390485 0.79390485 0.76871068 0.76871068 0.75689295
## [31] 0.75591538 0.75591538 0.75462884 0.75462884 0.75363509 0.75363509
## [37] 0.74120760 0.74120760 0.73545145 0.73545145 0.15559327 0.06386759(ammonia_resid_sd <- sd(resid(var_res_fit)[,1]))## [1] 1.616646(hydrogensulfide_resid_sd <- sd(resid(var_res_fit)[,2]))## [1] 16.47453(carbondioxide_resid_sd <- sd(resid(var_res_fit)[,3]))## [1] 13.21803\[ \begin{equation} \begin{cases} a_{t}\ =\ 0.23a_{t-1} + 0.004h_{t-1}\ + 0.01c_{t-1} + 0.13a_{t-2} - 0.005h_{t-2} +\\ \ \ \ \ \ \ \ \ \ +\ 0.15a_{t-3} - 0.003h_{t-2} + 0.1a_{t-4} + 0.08a_{t-5} - 0.004c_{t-5} +\\ \ \ \ \ \ \ \ \ \ +\ 0.06a_{t-6} + 0.03a_{t-8} + 0.003c_{t-8} + 0.05a_{t-9} + 0.06a_{t-10} - 0.006h_{t-10} +\\ \ \ \ \ \ \ \ \ \ +\ 0.01h_{t-11} + 0.08a_{t-12} + 0.003h_{t-13} - 0.004c_{t-13} + 0.03a_{t-14} +\\ \ \ \ \ \ \ \ \ \ +\ 0.01h_{t-14} - 0.00004t - 0.6(t==2:15am) + 0.56(t==5:30am) +\\ \ \ \ \ \ \ \ \ \ +\ 0.57(t==7:30am) + 0.49(t==3:00pm) - 0.46(t==7:15pm)\\ \ \ \ \ \ \ \ \ \ +\ 0.43(t==9:00pm) + 0.57(t==9:45am) - 0.5(t==11:15pm) + u_{t} \\ \\ u_{t}\ =\ \sigma \epsilon_{t}, \ \ \ \ \sigma\ =\ 1.63, \ \ \ \ \ \epsilon_{t}=N(0,1) \end{cases} \end{equation} \]
\[ \begin{equation} \begin{cases} h_{t}\ =\ 0.09h_{t-1} - 0.75a_{t-2} + 0.30a_{t-3} + 0.1h_{t-3} +0.03c_{t-3} - 0.05h_{t-4} + \\ \ \ \ \ \ \ \ \ \ \ + 0.03c_{t-4} - 0.05c_{t-5} - 0.02h_{t-6} + 0.38a_{t-7} + 0.17h_{t-7} + \\ \ \ \ \ \ \ \ \ \ \ - 0.05h{t-8} + 0.54a_{t-9} + 0.04a_{t-10} - 0.05h_{t-10} - 0.44a_{t-11} + \\ \ \ \ \ \ \ \ \ \ \ + 0.09h_{t-11} + 0.04h_{t-12} - 0.35a_{t-14} + 0.13h_{t-14} - 12.56 + \\ \ \ \ \ \ \ \ \ \ \ + 7.7(t==3:15am) + 4.64(t==9:30am) + 5.56(t==11:00am) + 4.96(t==8:45pm)\\ \ \ \ \ \ \ \ \ \ \ + 9.17(t==9:45pm) + 4.51(t==10pm) + 9.21(t==10:45pm) + 11.94(t==11:30pm)\\ \ \ \ \ \ \ \ \ \ \ + 5.19(t==11:45pm) + u_{t} \\ \\ u_{t}\ =\ \sigma \epsilon_{t}, \ \ \ \ \sigma\ =\ 16.63, \ \ \ \ \ \epsilon_{t}=N(0,1) \end{cases} \end{equation} \]
\[ \begin{equation} \begin{cases} c_{t}\ =\ -0.37a_{t-1} + 0.64c_{t-1} + 0.22c_{t-2} - 0.02h_{t-3} + 0.46a_{t-4} - 0.04h_{t-4}\\ \ \ \ \ \ \ \ \ \ \ \ \ + 0.02h_{t-5} + 0.04c_{t-5} + 0.02h{t-6} + 0.04c_{t-6} - 0.02h_{t-9}\\ \ \ \ \ \ \ \ \ \ \ \ \ - 0.02c_{t-9} - 0.03h_{t-11} + 0.02c_{t-14} + 30.04 - 0.00028t\\ \ \ \ \ \ \ \ \ \ \ \ \ - 4.17(t==2:00am) + 4.18(t==7:45am) + 5.32(t==3:30pm) \\ \ \ \ \ \ \ \ \ \ \ \ \ - 3.52(t==11:45pm) + u_{t} \\ \\ u_{t}\ =\ \sigma \epsilon_{t}, \ \ \ \ \sigma\ =\ 13.22, \ \ \ \ \ \epsilon_{t}=N(0,1) \end{cases} \end{equation} \]
var_res_stability <- stability(var_res_fit)
plot(var_res_stability$stability$ammonia)
sctest(var_res_stability$stability$ammonia)##
## OLS-based CUSUM test
##
## data: var_res_stability$stability$ammonia
## S0 = 1.3875, p-value = 0.04255plot(var_res_stability$stability$hydrogensulfide)
sctest(var_res_stability$stability$hydrogensulfide)##
## OLS-based CUSUM test
##
## data: var_res_stability$stability$hydrogensulfide
## S0 = 0.82081, p-value = 0.5107plot(var_res_stability$stability$carbondioxide)
sctest(var_res_stability$stability$carbondioxide)##
## OLS-based CUSUM test
##
## data: var_res_stability$stability$carbondioxide
## S0 = 1.113, p-value = 0.1678normality(var_res_fit, multivariate.only = TRUE)## $JB
##
## JB-Test (multivariate)
##
## data: Residuals of VAR object x
## Chi-squared = 34774000, df = 6, p-value < 2.2e-16
##
##
## $Skewness
##
## Skewness only (multivariate)
##
## data: Residuals of VAR object x
## Chi-squared = 271940, df = 3, p-value < 2.2e-16
##
##
## $Kurtosis
##
## Kurtosis only (multivariate)
##
## data: Residuals of VAR object x
## Chi-squared = 34502000, df = 3, p-value < 2.2e-16serial.test(var_res_fit, lags.pt = 14, lags.bg = 14, type = "PT.asymptotic")##
## Portmanteau Test (asymptotic)
##
## data: Residuals of VAR object var_res_fit
## Chi-squared = 123.31, df = 0, p-value < 2.2e-16serial.test(var_res_fit, lags.pt = 14, lags.bg = 14, type = "BG")##
## Breusch-Godfrey LM test
##
## data: Residuals of VAR object var_res_fit
## Chi-squared = 377.24, df = 126, p-value < 2.2e-16arch.test(var_res_fit, lags.multi = 14)##
## ARCH (multivariate)
##
## data: Residuals of VAR object var_res_fit
## Chi-squared = 4557, df = 504, p-value < 2.2e-16ammonia_irf <- irf(var_res_fit, impulse = "hydrogensulfide", response = c("ammonia"), boot = TRUE)
plot(ammonia_irf)
ammonia_irf <- irf(var_res_fit, impulse = "carbondioxide", response = c("ammonia"), boot = TRUE)
plot(ammonia_irf)
causality(var_res_fit, cause=c("hydrogensulfide", "carbondioxide"), vcov.=vcovHC(var_res_fit))## $Granger
##
## Granger causality H0: hydrogensulfide carbondioxide do not
## Granger-cause ammonia
##
## data: VAR object var_res_fit
## F-Test = 1.9831, df1 = 28, df2 = 17778, p-value = 0.001487
##
##
## $Instant
##
## H0: No instantaneous causality between: hydrogensulfide
## carbondioxide and ammonia
##
## data: VAR object var_res_fit
## Chi-squared = 188.7, df = 2, p-value < 2.2e-16causality(var_res_fit, cause="ammonia", vcov.=vcovHC(var_res_fit))## $Granger
##
## Granger causality H0: ammonia do not Granger-cause
## hydrogensulfide carbondioxide
##
## data: VAR object var_res_fit
## F-Test = 1.3269, df1 = 28, df2 = 17778, p-value = 0.1158
##
##
## $Instant
##
## H0: No instantaneous causality between: ammonia and
## hydrogensulfide carbondioxide
##
## data: VAR object var_res_fit
## Chi-squared = 188.7, df = 2, p-value < 2.2e-16causality(var_res_fit, cause="hydrogensulfide", vcov.=vcovHC(var_res_fit))## $Granger
##
## Granger causality H0: hydrogensulfide do not Granger-cause
## ammonia carbondioxide
##
## data: VAR object var_res_fit
## F-Test = 1.6494, df1 = 28, df2 = 17778, p-value = 0.01679
##
##
## $Instant
##
## H0: No instantaneous causality between: hydrogensulfide and
## ammonia carbondioxide
##
## data: VAR object var_res_fit
## Chi-squared = 54.585, df = 2, p-value = 1.403e-12causality(var_res_fit, cause="carbondioxide", vcov.=vcovHC(var_res_fit))## $Granger
##
## Granger causality H0: carbondioxide do not Granger-cause ammonia
## hydrogensulfide
##
## data: VAR object var_res_fit
## F-Test = 1.8924, df1 = 28, df2 = 17778, p-value = 0.002979
##
##
## $Instant
##
## H0: No instantaneous causality between: carbondioxide and ammonia
## hydrogensulfide
##
## data: VAR object var_res_fit
## Chi-squared = 139.61, df = 2, p-value < 2.2e-16plot(predict(var_res_fit, n.ahead=24, ci=0.9), xlim=c(5500,6100))
Utilizziamo i test di Phillips-Ouliaris e Johansen. La Null Hypothesis e’ di non co-integrazione.
po.test(gas[,c("ammonia", "carbondioxide")])##
## Phillips-Ouliaris Cointegration Test
##
## data: gas[, c("ammonia", "carbondioxide")]
## Phillips-Ouliaris demeaned = -4207.8, Truncation lag parameter =
## 60, p-value = 0.01po.test(gas)##
## Phillips-Ouliaris Cointegration Test
##
## data: gas
## Phillips-Ouliaris demeaned = -4484.6, Truncation lag parameter =
## 60, p-value = 0.01##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: trace statistic , without linear trend and constant in cointegration
##
## Eigenvalues (lambda):
## [1] 4.164497e-02 1.837614e-02 6.938894e-18
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 1 | 112.71 7.52 9.24 12.97
## r = 0 | 371.21 17.85 19.96 24.60
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## ammonia.l2 carbondioxide.l2 constant
## ammonia.l2 1.00000000 1.0000000 1.0000000
## carbondioxide.l2 -0.09189339 0.1813385 0.4087976
## constant 44.08172340 -108.8527985 -1507.7442444
##
## Weights W:
## (This is the loading matrix)
##
## ammonia.l2 carbondioxide.l2 constant
## ammonia.d -0.08617893 -0.0153477 2.124647e-19
## carbondioxide.d 0.33095254 -0.1794838 -1.322624e-18##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: trace statistic , without linear trend and constant in cointegration
##
## Eigenvalues (lambda):
## [1] 3.170234e-01 4.119962e-02 1.828980e-02 -2.574980e-19
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 2 | 112.18 7.52 9.24 12.97
## r <= 1 | 367.85 17.85 19.96 24.60
## r = 0 | 2684.98 32.00 34.91 41.07
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## ammonia.l2 hydrogensulfide.l2 carbondioxide.l2
## ammonia.l2 1.00000000 1.00000000 1.00000000
## hydrogensulfide.l2 -2.71467611 -0.01237143 -0.01394933
## carbondioxide.l2 0.05175219 -0.09221790 0.17389290
## constant -32.68892294 44.28000048 -104.66629133
## constant
## ammonia.l2 1.0000000
## hydrogensulfide.l2 -0.0188399
## carbondioxide.l2 0.3989170
## constant -1470.5043533
##
## Weights W:
## (This is the loading matrix)
##
## ammonia.l2 hydrogensulfide.l2 carbondioxide.l2
## ammonia.d 0.002117407 -0.08487263 -0.015962257
## hydrogensulfide.d 0.339872505 -0.07123760 -0.004885808
## carbondioxide.d 0.002634336 0.33845996 -0.183326908
## constant
## ammonia.d -1.453566e-19
## hydrogensulfide.d 8.520800e-20
## carbondioxide.d -2.026829e-18gas_vecm <- ca.jo(gas, ecdet = "trend", type = "trace", K = 14, spec="longrun", season = 96)
summary(gas_vecm)##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: trace statistic , with linear trend in cointegration
##
## Eigenvalues (lambda):
## [1] 3.634617e-02 1.355339e-02 4.033288e-03 -1.214826e-18
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 2 | 24.51 10.49 12.25 16.26
## r <= 1 | 107.27 22.76 25.32 30.45
## r = 0 | 331.82 39.06 42.44 48.45
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## ammonia.l14 hydrogensulfide.l14 carbondioxide.l14
## ammonia.l14 1.0000000000 1.000000000 1.000000000
## hydrogensulfide.l14 -1.9469000510 -0.056738909 -0.001212987
## carbondioxide.l14 0.0490428594 -0.318504170 -0.011630208
## trend.l14 0.0006566638 -0.001130951 0.001121571
## trend.l14
## ammonia.l14 1.00000000
## hydrogensulfide.l14 -0.02357937
## carbondioxide.l14 0.03009678
## trend.l14 0.02336693
##
## Weights W:
## (This is the loading matrix)
##
## ammonia.l14 hydrogensulfide.l14 carbondioxide.l14
## ammonia.d -0.0004875185 -0.003656281 -0.03571476
## hydrogensulfide.d 0.2829043125 -0.028782315 -0.03352622
## carbondioxide.d 0.0215571394 0.163204199 -0.09481076
## trend.l14
## ammonia.d -5.562107e-18
## hydrogensulfide.d -2.174549e-16
## carbondioxide.d 8.573699e-17(vec2var_res <- vec2var(gas_vecm))##
## Coefficient matrix of lagged endogenous variables:
##
## A1:
## ammonia.l1 hydrogensulfide.l1 carbondioxide.l1
## ammonia 0.2354099 0.003898711 0.009084324
## hydrogensulfide 0.1751574 0.095554399 0.004434282
## carbondioxide -0.4460881 0.017262244 0.657283091
##
##
## A2:
## ammonia.l2 hydrogensulfide.l2 carbondioxide.l2
## ammonia 0.1289521 -0.005355138 0.002368705
## hydrogensulfide -0.7905082 -0.020588705 -0.032579193
## carbondioxide 0.1046958 -0.020171187 0.224849717
##
##
## A3:
## ammonia.l3 hydrogensulfide.l3 carbondioxide.l3
## ammonia 0.1500075 0.001084417 -0.003860731
## hydrogensulfide 0.2243934 0.111567408 0.050330609
## carbondioxide -0.2380470 -0.027005754 0.020769700
##
##
## A4:
## ammonia.l4 hydrogensulfide.l4 carbondioxide.l4
## ammonia 0.10145525 -0.0003932661 -0.001124837
## hydrogensulfide 0.06291896 -0.0521061264 0.037520236
## carbondioxide 0.42509293 -0.0392249240 -0.004334919
##
##
## A5:
## ammonia.l5 hydrogensulfide.l5 carbondioxide.l5
## ammonia 0.08075696 -0.0007020275 -0.003663352
## hydrogensulfide 0.04909186 0.0036658361 -0.056774899
## carbondioxide 0.05074426 0.0261008354 0.035328042
##
##
## A6:
## ammonia.l6 hydrogensulfide.l6 carbondioxide.l6
## ammonia 0.06021504 0.001681929 -0.0005970547
## hydrogensulfide 0.01204024 -0.024816779 0.0075581850
## carbondioxide 0.19962173 0.020295124 0.0517951413
##
##
## A7:
## ammonia.l7 hydrogensulfide.l7 carbondioxide.l7
## ammonia 0.02640210 0.0002101925 -0.001032216
## hydrogensulfide 0.35558546 0.1800676420 -0.020816085
## carbondioxide 0.05906567 0.0059340968 -0.011009355
##
##
## A8:
## ammonia.l8 hydrogensulfide.l8 carbondioxide.l8
## ammonia 0.03288288 0.001114234 0.003324734
## hydrogensulfide 0.17757572 -0.054754440 0.020785008
## carbondioxide -0.20014528 0.001416897 0.003848144
##
##
## A9:
## ammonia.l9 hydrogensulfide.l9 carbondioxide.l9
## ammonia 0.05091685 0.0006874013 0.001551782
## hydrogensulfide 0.52566854 -0.0016941769 0.022651693
## carbondioxide 0.11691619 -0.0228880733 -0.030414059
##
##
## A10:
## ammonia.l10 hydrogensulfide.l10 carbondioxide.l10
## ammonia 0.06768477 -0.006159189 -0.0009078766
## hydrogensulfide 0.48652956 -0.051232457 -0.0372388048
## carbondioxide -0.12707533 0.016961274 0.0189429324
##
##
## A11:
## ammonia.l11 hydrogensulfide.l11 carbondioxide.l11
## ammonia -0.01071455 -0.01042746 0.000287801
## hydrogensulfide -0.46507724 0.09499211 -0.010825472
## carbondioxide 0.13631352 -0.03461657 0.001003019
##
##
## A12:
## ammonia.l12 hydrogensulfide.l12 carbondioxide.l12
## ammonia 0.08924602 0.002180133 -0.001321279
## hydrogensulfide 0.10451989 0.048424284 0.018588583
## carbondioxide -0.15325123 0.010522867 -0.018051252
##
##
## A13:
## ammonia.l13 hydrogensulfide.l13 carbondioxide.l13
## ammonia 0.012613807 0.002260957 -0.002947948
## hydrogensulfide -0.272196494 -0.016901116 0.004771959
## carbondioxide 0.003594361 0.003740416 0.006465014
##
##
## A14:
## ammonia.l14 hydrogensulfide.l14 carbondioxide.l14
## ammonia -0.02631619 0.0108682535 -0.001185960
## hydrogensulfide -0.36279478 0.1370357003 0.005468335
## carbondioxide 0.09011954 -0.0002968451 0.044582008
##
##
## Coefficient matrix of deterministic regressor(s).
##
## constant trend.l14 sd1 sd2
## ammonia 0.01385935 -3.201358e-07 -0.1522085 0.1804174
## hydrogensulfide -9.66906393 1.857730e-04 -1.0843149 -2.3932963
## carbondioxide -0.77575598 1.415579e-05 -0.5552372 -1.2806025
## sd3 sd4 sd5 sd6 sd7
## ammonia -0.16250413 0.1130428 0.1024962 0.4697432 0.1464656
## hydrogensulfide -0.44713960 -0.3302213 2.9097599 -1.6728565 -0.1498103
## carbondioxide 0.02288145 0.7164059 1.8327884 0.5508000 -1.9316969
## sd8 sd9 sd10 sd11 sd12
## ammonia 0.21052922 -0.4761969 0.3839420 0.2729766 -0.00710544
## hydrogensulfide 0.09212816 2.6782872 -1.4097127 -0.4053443 -2.73423482
## carbondioxide -4.79436811 -2.8809862 0.5606761 1.0515582 -2.60638765
## sd13 sd14 sd15 sd16 sd17
## ammonia -0.02823491 0.1244017 -0.2021090 0.2317202 0.1131266
## hydrogensulfide 7.79644796 -0.6867082 0.9639996 0.8108020 1.5947227
## carbondioxide 0.46031098 1.5340869 0.2023734 1.0119119 -0.4900680
## sd18 sd19 sd20 sd21 sd22
## ammonia 0.05753545 0.21226877 -0.09304993 0.1768923 0.7149108
## hydrogensulfide 2.06290697 0.06783507 -1.33794915 0.3610656 0.1742328
## carbondioxide 0.02535349 1.63357943 -1.19139084 -2.7847439 -2.4563807
## sd23 sd24 sd25 sd26 sd27
## ammonia -0.1079527 0.3056529 -0.0635461 0.30836587 0.2059225
## hydrogensulfide 0.9026026 1.0141477 -0.3892589 0.62852201 1.4853055
## carbondioxide -2.4781948 2.1109744 -0.2166917 0.01527936 -0.2897014
## sd28 sd29 sd30 sd31 sd32
## ammonia 0.0264630 0.2020182 0.71478661 0.1706784 0.07948712
## hydrogensulfide -0.5287123 0.5901102 -0.14992291 -0.3983468 0.55778985
## carbondioxide -0.1003034 0.8485584 0.06944746 4.0071796 -0.76091848
## sd33 sd34 sd35 sd36 sd37
## ammonia 0.05934073 -0.1099529 0.1028578 -0.01394281 0.2864333
## hydrogensulfide 1.10065886 -0.9355637 0.9826651 -0.37523058 -0.5082938
## carbondioxide 0.27283318 0.4828906 -3.1095933 -3.34180817 0.7216260
## sd38 sd39 sd40 sd41 sd42
## ammonia 0.1400465 -0.02143027 0.2671349 0.2467514 0.11617982
## hydrogensulfide 4.5098332 -0.78884219 1.0931525 -1.1870181 1.71922509
## carbondioxide 0.1989907 -2.30742488 -1.0773812 -1.5371084 0.02729479
## sd43 sd44 sd45 sd46 sd47
## ammonia 0.02215667 0.217772 0.1315505 0.1189649 0.3222925
## hydrogensulfide -0.02629371 5.833406 -1.0745365 0.4390557 -1.2308737
## carbondioxide 0.29094630 1.467551 -0.3148527 -1.3146695 -1.3245417
## sd48 sd49 sd50 sd51 sd52
## ammonia 0.1551531 0.2531826 0.09828939 0.1735216 0.2667105
## hydrogensulfide 0.4624296 -0.8014935 0.70200106 -0.9421571 -0.2802047
## carbondioxide 0.7354533 1.2210218 -0.75560167 1.7450999 0.5258438
## sd53 sd54 sd55 sd56 sd57
## ammonia 0.1371479 0.1724188 -0.03580989 -0.04238372 0.2389521
## hydrogensulfide -0.1771672 0.2238351 0.29303149 -0.39399982 -0.1786948
## carbondioxide -3.0978824 -1.0490150 2.21778479 -2.08662578 0.6670631
## sd58 sd59 sd60 sd61 sd62
## ammonia 0.5426416 0.1905049 0.6209294 -0.2304553 0.2710306
## hydrogensulfide 2.4935452 -0.5011652 1.3406474 -0.5853145 0.5972393
## carbondioxide 0.1981159 -1.5445210 0.3143734 -0.3027751 5.2278729
## sd63 sd64 sd65 sd66 sd67
## ammonia 0.2555313 0.1989378 0.19105883 -0.2336031 0.1580402
## hydrogensulfide -0.3147949 0.6878334 0.04343838 0.6516539 -0.6417860
## carbondioxide 1.1659897 -1.6330344 -0.41619471 -1.5577549 0.4567582
## sd68 sd69 sd70 sd71 sd72
## ammonia 0.03760525 0.03936553 0.09005637 0.17513009 0.2548263
## hydrogensulfide 1.00485733 -0.70437520 -0.12104266 -0.38476134 -0.3937060
## carbondioxide 0.70761643 -0.53732119 2.66920972 0.02518144 1.1933465
## sd73 sd74 sd75 sd76 sd77
## ammonia 0.24010483 0.1285136 -0.0869835 0.00285022 -0.3280251
## hydrogensulfide -0.05474672 0.8128941 0.1095724 0.40196296 -0.1586556
## carbondioxide -0.42565600 -0.7791407 -3.0294100 -0.56904948 -0.3347144
## sd78 sd79 sd80 sd81 sd82
## ammonia -0.09037394 0.3034828 0.1297794 -0.1710675 0.2593676
## hydrogensulfide 1.68520179 -0.6781066 -0.6420060 -0.1836844 -0.2035127
## carbondioxide -0.08089691 -0.1220891 -1.0711430 -1.1311780 0.5686033
## sd83 sd84 sd85 sd86 sd87
## ammonia 0.3021801 0.5875631 0.3212924 0.7015305 0.3708868
## hydrogensulfide 5.0797315 0.1436071 0.1448829 2.1623581 9.2616094
## carbondioxide -0.4456442 -0.2054663 2.7243675 0.6481400 1.0171343
## sd88 sd89 sd90 sd91 sd92
## ammonia 0.1324237 0.4858997 0.2804017 0.6107419 -0.1584471
## hydrogensulfide 4.5200609 1.2728897 -2.0535643 9.1134096 0.1134191
## carbondioxide -3.5001239 -0.3045013 0.7821833 -0.4114435 -1.8764777
## sd93 sd94 sd95
## ammonia -0.3928071 0.2110355 0.08204357
## hydrogensulfide 1.1149790 11.3611678 4.99473765
## carbondioxide -3.9500169 0.5776481 -0.75472101fit_ts <- ts(fitted(vec2var_res)[,1], frequency = frequency(ammonia_ts))
rmse <- accuracy(ammonia_ts, fit_ts)[,"RMSE"]
plot(ammonia_ts, ylab = 'ammonia', lwd = 1, lty ='dashed'); grid(); lines(fit_ts, col = 'red', lty = 'dashed')
legend("topright", c("ammonia", "fitted"), fill=c("black", "red")); rmse <- accuracy(ammonia_ts, fit_ts)[,"RMSE"]
mes <- paste("RMSE = ", round(rmse, 2), sep=''); text(60, 30, mes);
fit_ts <- ts(fitted(vec2var_res)[,2], frequency = frequency(hydrogensulfide_ts))
rmse <- accuracy(hydrogensulfide_ts, fit_ts)[,"RMSE"]
plot(hydrogensulfide_ts, ylab = 'hydrogensulfide', lwd = 1, lty ='dashed'); grid(); lines(fit_ts, col = 'red', lty = 'dashed')
legend("topright", c("hydrogensulfide", "fitted"), fill=c("black", "red")); rmse <- accuracy(hydrogensulfide_ts, fit_ts)[,"RMSE"]
mes <- paste("RMSE = ", round(rmse, 2), sep=''); text(60, 200, mes)
fit_ts <- ts(fitted(vec2var_res)[,3], frequency = frequency(carbondioxide_ts))
rmse <- accuracy(carbondioxide_ts, fit_ts)[,"RMSE"]
plot(carbondioxide_ts, ylab = 'carbondioxide', lwd = 1, lty ='dashed'); grid(); lines(fit_ts, col = 'red', lty = 'dashed')
legend("topright", c("carbondioxide", "fitted"), fill=c("black", "red")); rmse <- accuracy(carbondioxide_ts, fit_ts)[,"RMSE"]
mes <- paste("RMSE = ", round(rmse, 2), sep=''); text(60, 700, mes)
Viene definito un modello predittivo basato su algoritmo di regressione logistica.
ammonia_outlier <- ammonia_outlier_absolute
ammonia_outliers_indication <- rep(0, length(ammonia_ts))
ammonia_outliers_indication[ammonia_outlier$iRight] <- 1
ammonia_outliers_df <- data.frame(ammonia_outliers_indication,
Hmisc::Lag(ammonia_ts, shift = 4),
Hmisc::Lag(carbondioxide_ts, shift = 4),
Hmisc::Lag(ammonia_ts, shift = 5),
Hmisc::Lag(carbondioxide_ts, shift = 5)
)
ammonia_outliers_df <- ammonia_outliers_df[complete.cases(ammonia_outliers_df),]
colnames(ammonia_outliers_df) <- c("ammonia_peak",
"ammonia_4",
"carbondioxide_4",
"ammonia_5",
"carbondioxide_5"
)
ammonia_outliers_df$ammonia_peak <- factor(ammonia_outliers_df$ammonia_peak)
levels(ammonia_outliers_df$ammonia_peak) <- c("F", "T")l <- nrow(ammonia_outliers_df)
nn <- 8
grp_set <- seq(1:(l/nn))
l_grp <- round(0.6 * length(grp_set))
grp_set_train <- sample(grp_set, l_grp, replace=FALSE)
train_idx <- do.call("c", lapply(grp_set_train, function(x) {((x-1)*nn):(x*nn)}))
trControl <- trainControl(method = "repeatedcv",
number = 10,
repeats = 5,
verboseIter = FALSE,
sampling = "smote",
classProbs = TRUE,
summaryFunction = twoClassSummary,
savePredictions = TRUE,
returnResamp = "all")glm_fit <- train(
ammonia_peak ~ .,
data = ammonia_outliers_df[train_idx,],
method = "glm",
trControl = trControl,
metric = 'ROC',
preProcess = c("BoxCox", "center", "scale")
)glm_fit## Generalized Linear Model
##
## 4095 samples
## 4 predictor
## 2 classes: 'F', 'T'
##
## Pre-processing: Box-Cox transformation (2), centered (4), scaled (4)
## Resampling: Cross-Validated (10 fold, repeated 5 times)
## Summary of sample sizes: 3685, 3686, 3685, 3685, 3686, 3686, ...
## Addtional sampling using SMOTE prior to pre-processing
##
## Resampling results:
##
## ROC Sens Spec
## 0.9740295 0.9591807 0.9307692summary(glm_fit$finalModel)##
## Call:
## NULL
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.1448 -0.1814 -0.0948 0.2280 3.4165
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.7992 0.1239 -6.453 1.10e-10 ***
## ammonia_4 1.0821 0.1655 6.538 6.23e-11 ***
## carbondioxide_4 1.1138 0.2707 4.114 3.89e-05 ***
## ammonia_5 2.8729 0.2044 14.054 < 2e-16 ***
## carbondioxide_5 -0.8757 0.2670 -3.280 0.00104 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2428.42 on 1777 degrees of freedom
## Residual deviance: 586.55 on 1773 degrees of freedom
## AIC: 596.55
##
## Number of Fisher Scoring iterations: 7varImp(glm_fit)## glm variable importance
##
## Overall
## ammonia_5 100.000
## ammonia_4 30.244
## carbondioxide_4 7.746
## carbondioxide_5 0.000train_set <- ammonia_outliers_df[train_idx,]
pred <- ROCR::prediction(as.vector(predict(glm_fit, train_set, type="prob")[,2]), train_set$ammonia_peak)
perf <- performance(pred, "tpr", "fpr")
plot(perf, colorize=TRUE, print.cutoffs.at=c(0.3, 0.5, 0.7))
glm_train_pred <- predict(glm_fit, train_set, type="prob")
pred_values <- ifelse(glm_train_pred[,2] > 0.8, "T", "F")
(cm_mat <- confusionMatrix(factor(pred_values), train_set$ammonia_peak))## Confusion Matrix and Statistics
##
## Reference
## Prediction F T
## F 3801 34
## T 40 220
##
## Accuracy : 0.9819
## 95% CI : (0.9774, 0.9858)
## No Information Rate : 0.938
## P-Value [Acc > NIR] : <2e-16
##
## Kappa : 0.8464
## Mcnemar's Test P-Value : 0.5611
##
## Sensitivity : 0.9896
## Specificity : 0.8661
## Pos Pred Value : 0.9911
## Neg Pred Value : 0.8462
## Prevalence : 0.9380
## Detection Rate : 0.9282
## Detection Prevalence : 0.9365
## Balanced Accuracy : 0.9279
##
## 'Positive' Class : F
## test_set <- ammonia_outliers_df[-train_idx,]
glm_test_pred <- predict(glm_fit, test_set, type="prob")
pred_values <- ifelse(glm_test_pred[,2] > 0.8, "T", "F")
(cm_mat <- confusionMatrix(factor(pred_values), test_set$ammonia_peak))## Confusion Matrix and Statistics
##
## Reference
## Prediction F T
## F 2100 30
## T 21 101
##
## Accuracy : 0.9774
## 95% CI : (0.9703, 0.9831)
## No Information Rate : 0.9418
## P-Value [Acc > NIR] : 2.861e-16
##
## Kappa : 0.7864
## Mcnemar's Test P-Value : 0.2626
##
## Sensitivity : 0.9901
## Specificity : 0.7710
## Pos Pred Value : 0.9859
## Neg Pred Value : 0.8279
## Prevalence : 0.9418
## Detection Rate : 0.9325
## Detection Prevalence : 0.9458
## Balanced Accuracy : 0.8805
##
## 'Positive' Class : F
## l <- nrow(gas); s <- l/4; ll <- as.list(1:s); ll2 <- lapply(ll, function(x){ rep(x, 4)}); sample <- do.call(c, ll2); ls <- length(sample)
df <- data.frame(gas$ammonia[1:ls], sample); colnames(df) <- c("ammonia", "sample")
ammonia_g = with(df, qcc.groups(ammonia, sample))
ammonia_qchart <- qcc(ammonia_g[961:985,], type="xbar", newdata=ammonia_g[986:1034], nsigmas=3)
Riferimento: http://www.statsoft.com/textbook/quality-control-charts