R version 2.13.0 (2011-04-13)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-pc-linux-gnu (32-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(12008,9169,8788,8417,8247,8197,8236,8253,7733,8366,8626,8863,10102,8463,9114,8563,8872,8301,8301,8278,7736,7973,8268,9476,11100,8962,9173,8738,8459,8078,8411,8291,7810,8616,8312,9692,9911,8915,9452,9112,8472,8230,8384,8625,8221,8649,8625,10443,10357,8586,8892,8329,8101,7922,8120,7838,7735,8406,8209,9451,10041,9411,10405,8467,8464,8102,7627,7513,7510,8291,8064,9383,9706,8579,9474,8318,8213,8059,9111,7708,7680,8014,8007,8718,9486,9113,9025,8476,7952,7759,7835,7600,7651,8319,8812,8630),dim=c(1,96),dimnames=list(c('Maandelijkse_sterftes'),1:96))
> y <- array(NA,dim=c(1,96),dimnames=list(c('Maandelijkse_sterftes'),1:96))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
Maandelijkse_sterftes M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 12008 1 0 0 0 0 0 0 0 0 0 0 1
2 9169 0 1 0 0 0 0 0 0 0 0 0 2
3 8788 0 0 1 0 0 0 0 0 0 0 0 3
4 8417 0 0 0 1 0 0 0 0 0 0 0 4
5 8247 0 0 0 0 1 0 0 0 0 0 0 5
6 8197 0 0 0 0 0 1 0 0 0 0 0 6
7 8236 0 0 0 0 0 0 1 0 0 0 0 7
8 8253 0 0 0 0 0 0 0 1 0 0 0 8
9 7733 0 0 0 0 0 0 0 0 1 0 0 9
10 8366 0 0 0 0 0 0 0 0 0 1 0 10
11 8626 0 0 0 0 0 0 0 0 0 0 1 11
12 8863 0 0 0 0 0 0 0 0 0 0 0 12
13 10102 1 0 0 0 0 0 0 0 0 0 0 13
14 8463 0 1 0 0 0 0 0 0 0 0 0 14
15 9114 0 0 1 0 0 0 0 0 0 0 0 15
16 8563 0 0 0 1 0 0 0 0 0 0 0 16
17 8872 0 0 0 0 1 0 0 0 0 0 0 17
18 8301 0 0 0 0 0 1 0 0 0 0 0 18
19 8301 0 0 0 0 0 0 1 0 0 0 0 19
20 8278 0 0 0 0 0 0 0 1 0 0 0 20
21 7736 0 0 0 0 0 0 0 0 1 0 0 21
22 7973 0 0 0 0 0 0 0 0 0 1 0 22
23 8268 0 0 0 0 0 0 0 0 0 0 1 23
24 9476 0 0 0 0 0 0 0 0 0 0 0 24
25 11100 1 0 0 0 0 0 0 0 0 0 0 25
26 8962 0 1 0 0 0 0 0 0 0 0 0 26
27 9173 0 0 1 0 0 0 0 0 0 0 0 27
28 8738 0 0 0 1 0 0 0 0 0 0 0 28
29 8459 0 0 0 0 1 0 0 0 0 0 0 29
30 8078 0 0 0 0 0 1 0 0 0 0 0 30
31 8411 0 0 0 0 0 0 1 0 0 0 0 31
32 8291 0 0 0 0 0 0 0 1 0 0 0 32
33 7810 0 0 0 0 0 0 0 0 1 0 0 33
34 8616 0 0 0 0 0 0 0 0 0 1 0 34
35 8312 0 0 0 0 0 0 0 0 0 0 1 35
36 9692 0 0 0 0 0 0 0 0 0 0 0 36
37 9911 1 0 0 0 0 0 0 0 0 0 0 37
38 8915 0 1 0 0 0 0 0 0 0 0 0 38
39 9452 0 0 1 0 0 0 0 0 0 0 0 39
40 9112 0 0 0 1 0 0 0 0 0 0 0 40
41 8472 0 0 0 0 1 0 0 0 0 0 0 41
42 8230 0 0 0 0 0 1 0 0 0 0 0 42
43 8384 0 0 0 0 0 0 1 0 0 0 0 43
44 8625 0 0 0 0 0 0 0 1 0 0 0 44
45 8221 0 0 0 0 0 0 0 0 1 0 0 45
46 8649 0 0 0 0 0 0 0 0 0 1 0 46
47 8625 0 0 0 0 0 0 0 0 0 0 1 47
48 10443 0 0 0 0 0 0 0 0 0 0 0 48
49 10357 1 0 0 0 0 0 0 0 0 0 0 49
50 8586 0 1 0 0 0 0 0 0 0 0 0 50
51 8892 0 0 1 0 0 0 0 0 0 0 0 51
52 8329 0 0 0 1 0 0 0 0 0 0 0 52
53 8101 0 0 0 0 1 0 0 0 0 0 0 53
54 7922 0 0 0 0 0 1 0 0 0 0 0 54
55 8120 0 0 0 0 0 0 1 0 0 0 0 55
56 7838 0 0 0 0 0 0 0 1 0 0 0 56
57 7735 0 0 0 0 0 0 0 0 1 0 0 57
58 8406 0 0 0 0 0 0 0 0 0 1 0 58
59 8209 0 0 0 0 0 0 0 0 0 0 1 59
60 9451 0 0 0 0 0 0 0 0 0 0 0 60
61 10041 1 0 0 0 0 0 0 0 0 0 0 61
62 9411 0 1 0 0 0 0 0 0 0 0 0 62
63 10405 0 0 1 0 0 0 0 0 0 0 0 63
64 8467 0 0 0 1 0 0 0 0 0 0 0 64
65 8464 0 0 0 0 1 0 0 0 0 0 0 65
66 8102 0 0 0 0 0 1 0 0 0 0 0 66
67 7627 0 0 0 0 0 0 1 0 0 0 0 67
68 7513 0 0 0 0 0 0 0 1 0 0 0 68
69 7510 0 0 0 0 0 0 0 0 1 0 0 69
70 8291 0 0 0 0 0 0 0 0 0 1 0 70
71 8064 0 0 0 0 0 0 0 0 0 0 1 71
72 9383 0 0 0 0 0 0 0 0 0 0 0 72
73 9706 1 0 0 0 0 0 0 0 0 0 0 73
74 8579 0 1 0 0 0 0 0 0 0 0 0 74
75 9474 0 0 1 0 0 0 0 0 0 0 0 75
76 8318 0 0 0 1 0 0 0 0 0 0 0 76
77 8213 0 0 0 0 1 0 0 0 0 0 0 77
78 8059 0 0 0 0 0 1 0 0 0 0 0 78
79 9111 0 0 0 0 0 0 1 0 0 0 0 79
80 7708 0 0 0 0 0 0 0 1 0 0 0 80
81 7680 0 0 0 0 0 0 0 0 1 0 0 81
82 8014 0 0 0 0 0 0 0 0 0 1 0 82
83 8007 0 0 0 0 0 0 0 0 0 0 1 83
84 8718 0 0 0 0 0 0 0 0 0 0 0 84
85 9486 1 0 0 0 0 0 0 0 0 0 0 85
86 9113 0 1 0 0 0 0 0 0 0 0 0 86
87 9025 0 0 1 0 0 0 0 0 0 0 0 87
88 8476 0 0 0 1 0 0 0 0 0 0 0 88
89 7952 0 0 0 0 1 0 0 0 0 0 0 89
90 7759 0 0 0 0 0 1 0 0 0 0 0 90
91 7835 0 0 0 0 0 0 1 0 0 0 0 91
92 7600 0 0 0 0 0 0 0 1 0 0 0 92
93 7651 0 0 0 0 0 0 0 0 1 0 0 93
94 8319 0 0 0 0 0 0 0 0 0 1 0 94
95 8812 0 0 0 0 0 0 0 0 0 0 1 95
96 8630 0 0 0 0 0 0 0 0 0 0 0 96
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) M1 M2 M3 M4 M5
9560.571 960.314 -474.578 -79.720 -813.362 -1014.130
M6 M7 M8 M9 M10 M11
-1276.397 -1100.039 -1335.681 -1585.198 -1011.216 -970.858
t
-4.233
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-680.15 -218.57 -19.82 124.50 1491.35
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9560.571 164.960 57.957 < 2e-16 ***
M1 960.314 203.618 4.716 9.59e-06 ***
M2 -474.578 203.501 -2.332 0.022120 *
M3 -79.720 203.395 -0.392 0.696101
M4 -813.362 203.300 -4.001 0.000136 ***
M5 -1014.130 203.216 -4.990 3.27e-06 ***
M6 -1276.397 203.144 -6.283 1.46e-08 ***
M7 -1100.039 203.082 -5.417 5.80e-07 ***
M8 -1335.681 203.032 -6.579 3.99e-09 ***
M9 -1585.198 202.993 -7.809 1.56e-11 ***
M10 -1011.216 202.965 -4.982 3.38e-06 ***
M11 -970.858 202.948 -4.784 7.38e-06 ***
t -4.233 1.507 -2.809 0.006187 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 405.9 on 83 degrees of freedom
Multiple R-squared: 0.7757, Adjusted R-squared: 0.7433
F-statistic: 23.93 on 12 and 83 DF, p-value: < 2.2e-16
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 0.9870200 0.02595997 0.012979986
[2,] 0.9932155 0.01356905 0.006784524
[3,] 0.9871334 0.02573330 0.012866648
[4,] 0.9767903 0.04641946 0.023209728
[5,] 0.9594030 0.08119398 0.040596990
[6,] 0.9354202 0.12915959 0.064579796
[7,] 0.9201474 0.15970527 0.079852635
[8,] 0.8884826 0.22303482 0.111517411
[9,] 0.9009656 0.19806887 0.099034435
[10,] 0.9058265 0.18834699 0.094173493
[11,] 0.8759878 0.24802450 0.124012248
[12,] 0.8605860 0.27882794 0.139413968
[13,] 0.8228778 0.35424443 0.177122217
[14,] 0.7674498 0.46510047 0.232550234
[15,] 0.7124055 0.57518896 0.287594482
[16,] 0.6500217 0.69995658 0.349978292
[17,] 0.5786927 0.84261470 0.421307348
[18,] 0.5126374 0.97472523 0.487362615
[19,] 0.4848100 0.96961997 0.515190015
[20,] 0.4323147 0.86462943 0.567685286
[21,] 0.4170191 0.83403821 0.582980897
[22,] 0.6568587 0.68628253 0.343141266
[23,] 0.6032830 0.79343397 0.396716984
[24,] 0.5892885 0.82142298 0.410711490
[25,] 0.5977951 0.80440971 0.402204853
[26,] 0.5301900 0.93962008 0.469810040
[27,] 0.4608798 0.92175965 0.539120173
[28,] 0.3941805 0.78836095 0.605819524
[29,] 0.4108589 0.82171771 0.589141143
[30,] 0.3897491 0.77949815 0.610250927
[31,] 0.3400590 0.68011798 0.659941008
[32,] 0.2847275 0.56945508 0.715272460
[33,] 0.6277139 0.74457222 0.372286108
[34,] 0.6649383 0.67012346 0.335061729
[35,] 0.6756024 0.64879515 0.324397577
[36,] 0.7619573 0.47608533 0.238042663
[37,] 0.7414211 0.51715778 0.258578888
[38,] 0.7229456 0.55410874 0.277054372
[39,] 0.6844361 0.63112777 0.315563884
[40,] 0.6420587 0.71588262 0.357941310
[41,] 0.6094836 0.78103285 0.390516424
[42,] 0.5436814 0.91263729 0.456318645
[43,] 0.4730491 0.94609821 0.526950893
[44,] 0.4315821 0.86316416 0.568417920
[45,] 0.3842568 0.76851353 0.615743237
[46,] 0.3750225 0.75004504 0.624977481
[47,] 0.3971447 0.79428932 0.602855340
[48,] 0.7663056 0.46738885 0.233694426
[49,] 0.7052513 0.58949734 0.294748672
[50,] 0.6657682 0.66846361 0.334231807
[51,] 0.6004337 0.79913261 0.399566305
[52,] 0.7542784 0.49144324 0.245721619
[53,] 0.7311383 0.53772333 0.268861665
[54,] 0.6811621 0.63767578 0.318837892
[55,] 0.5979370 0.80412606 0.402063032
[56,] 0.6017243 0.79655139 0.398275696
[57,] 0.6045651 0.79086984 0.395434920
[58,] 0.5534951 0.89300985 0.446504924
[59,] 0.5612802 0.87743969 0.438719845
[60,] 0.4929242 0.98584833 0.507075837
[61,] 0.4019149 0.80382980 0.598085101
[62,] 0.2982534 0.59650687 0.701746566
[63,] 0.2083213 0.41664252 0.791678739
[64,] 0.7786636 0.44267286 0.221336432
[65,] 0.6886518 0.62269642 0.311348211
> postscript(file="/var/wessaorg/rcomp/tmp/1gjku1322577850.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/2lj0a1322577850.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/3f3a51322577850.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/41gbn1322577850.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/59wrz1322577850.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> qqnorm(mysum$resid, main='Residual Normal Q-Q Plot')
> qqline(mysum$resid)
> grid()
> dev.off()
null device
1
> (myerror <- as.ts(mysum$resid))
Time Series:
Start = 1
End = 96
Frequency = 1
1 2 3 4 5 6
1491.3472222 91.4722222 -680.1527778 -313.2777778 -278.2777778 -61.7777778
7 8 9 10 11 12
-194.9027778 61.9722222 -204.2777778 -141.0277778 82.8472222 -646.7777778
13 14 15 16 17 18
-363.8591270 -563.7341270 -303.3591270 -116.4841270 397.5158730 93.0158730
19 20 21 22 23 24
-79.1091270 137.7658730 -150.4841270 -483.2341270 -224.3591270 17.0158730
25 26 27 28 29 30
684.9345238 -13.9404762 -193.5654762 109.3095238 35.3095238 -79.1904762
31 32 33 34 35 36
81.6845238 201.5595238 -25.6904762 210.5595238 -129.5654762 283.8095238
37 38 39 40 41 42
-453.2718254 -10.1468254 136.2281746 534.1031746 99.1031746 123.6031746
43 44 45 46 47 48
105.4781746 586.3531746 436.1031746 294.3531746 234.2281746 1085.6031746
49 50 51 52 53 54
43.5218254 -288.3531746 -372.9781746 -198.1031746 -221.1031746 -133.6031746
55 56 57 58 59 60
-107.7281746 -149.8531746 0.8968254 102.1468254 -130.9781746 144.3968254
61 62 63 64 65 66
-221.6845238 587.4404762 1190.8154762 -9.3095238 192.6904762 97.1904762
67 68 69 70 71 72
-549.9345238 -424.0595238 -173.3095238 37.9404762 -225.1845238 127.1904762
73 74 75 76 77 78
-505.8908730 -193.7658730 310.6091270 -107.5158730 -7.5158730 104.9841270
79 80 81 82 83 84
984.8591270 -178.2658730 47.4841270 -188.2658730 -231.3908730 -487.0158730
85 86 87 88 89 90
-675.0972222 391.0277778 -87.5972222 101.2777778 -217.7222222 -144.2222222
91 92 93 94 95 96
-240.3472222 -235.4722222 69.2777778 167.5277778 624.4027778 -524.2222222
> postscript(file="/var/wessaorg/rcomp/tmp/678n71322577850.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 96
Frequency = 1
lag(myerror, k = 1) myerror
0 1491.3472222 NA
1 91.4722222 1491.3472222
2 -680.1527778 91.4722222
3 -313.2777778 -680.1527778
4 -278.2777778 -313.2777778
5 -61.7777778 -278.2777778
6 -194.9027778 -61.7777778
7 61.9722222 -194.9027778
8 -204.2777778 61.9722222
9 -141.0277778 -204.2777778
10 82.8472222 -141.0277778
11 -646.7777778 82.8472222
12 -363.8591270 -646.7777778
13 -563.7341270 -363.8591270
14 -303.3591270 -563.7341270
15 -116.4841270 -303.3591270
16 397.5158730 -116.4841270
17 93.0158730 397.5158730
18 -79.1091270 93.0158730
19 137.7658730 -79.1091270
20 -150.4841270 137.7658730
21 -483.2341270 -150.4841270
22 -224.3591270 -483.2341270
23 17.0158730 -224.3591270
24 684.9345238 17.0158730
25 -13.9404762 684.9345238
26 -193.5654762 -13.9404762
27 109.3095238 -193.5654762
28 35.3095238 109.3095238
29 -79.1904762 35.3095238
30 81.6845238 -79.1904762
31 201.5595238 81.6845238
32 -25.6904762 201.5595238
33 210.5595238 -25.6904762
34 -129.5654762 210.5595238
35 283.8095238 -129.5654762
36 -453.2718254 283.8095238
37 -10.1468254 -453.2718254
38 136.2281746 -10.1468254
39 534.1031746 136.2281746
40 99.1031746 534.1031746
41 123.6031746 99.1031746
42 105.4781746 123.6031746
43 586.3531746 105.4781746
44 436.1031746 586.3531746
45 294.3531746 436.1031746
46 234.2281746 294.3531746
47 1085.6031746 234.2281746
48 43.5218254 1085.6031746
49 -288.3531746 43.5218254
50 -372.9781746 -288.3531746
51 -198.1031746 -372.9781746
52 -221.1031746 -198.1031746
53 -133.6031746 -221.1031746
54 -107.7281746 -133.6031746
55 -149.8531746 -107.7281746
56 0.8968254 -149.8531746
57 102.1468254 0.8968254
58 -130.9781746 102.1468254
59 144.3968254 -130.9781746
60 -221.6845238 144.3968254
61 587.4404762 -221.6845238
62 1190.8154762 587.4404762
63 -9.3095238 1190.8154762
64 192.6904762 -9.3095238
65 97.1904762 192.6904762
66 -549.9345238 97.1904762
67 -424.0595238 -549.9345238
68 -173.3095238 -424.0595238
69 37.9404762 -173.3095238
70 -225.1845238 37.9404762
71 127.1904762 -225.1845238
72 -505.8908730 127.1904762
73 -193.7658730 -505.8908730
74 310.6091270 -193.7658730
75 -107.5158730 310.6091270
76 -7.5158730 -107.5158730
77 104.9841270 -7.5158730
78 984.8591270 104.9841270
79 -178.2658730 984.8591270
80 47.4841270 -178.2658730
81 -188.2658730 47.4841270
82 -231.3908730 -188.2658730
83 -487.0158730 -231.3908730
84 -675.0972222 -487.0158730
85 391.0277778 -675.0972222
86 -87.5972222 391.0277778
87 101.2777778 -87.5972222
88 -217.7222222 101.2777778
89 -144.2222222 -217.7222222
90 -240.3472222 -144.2222222
91 -235.4722222 -240.3472222
92 69.2777778 -235.4722222
93 167.5277778 69.2777778
94 624.4027778 167.5277778
95 -524.2222222 624.4027778
96 NA -524.2222222
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 91.4722222 1491.3472222
[2,] -680.1527778 91.4722222
[3,] -313.2777778 -680.1527778
[4,] -278.2777778 -313.2777778
[5,] -61.7777778 -278.2777778
[6,] -194.9027778 -61.7777778
[7,] 61.9722222 -194.9027778
[8,] -204.2777778 61.9722222
[9,] -141.0277778 -204.2777778
[10,] 82.8472222 -141.0277778
[11,] -646.7777778 82.8472222
[12,] -363.8591270 -646.7777778
[13,] -563.7341270 -363.8591270
[14,] -303.3591270 -563.7341270
[15,] -116.4841270 -303.3591270
[16,] 397.5158730 -116.4841270
[17,] 93.0158730 397.5158730
[18,] -79.1091270 93.0158730
[19,] 137.7658730 -79.1091270
[20,] -150.4841270 137.7658730
[21,] -483.2341270 -150.4841270
[22,] -224.3591270 -483.2341270
[23,] 17.0158730 -224.3591270
[24,] 684.9345238 17.0158730
[25,] -13.9404762 684.9345238
[26,] -193.5654762 -13.9404762
[27,] 109.3095238 -193.5654762
[28,] 35.3095238 109.3095238
[29,] -79.1904762 35.3095238
[30,] 81.6845238 -79.1904762
[31,] 201.5595238 81.6845238
[32,] -25.6904762 201.5595238
[33,] 210.5595238 -25.6904762
[34,] -129.5654762 210.5595238
[35,] 283.8095238 -129.5654762
[36,] -453.2718254 283.8095238
[37,] -10.1468254 -453.2718254
[38,] 136.2281746 -10.1468254
[39,] 534.1031746 136.2281746
[40,] 99.1031746 534.1031746
[41,] 123.6031746 99.1031746
[42,] 105.4781746 123.6031746
[43,] 586.3531746 105.4781746
[44,] 436.1031746 586.3531746
[45,] 294.3531746 436.1031746
[46,] 234.2281746 294.3531746
[47,] 1085.6031746 234.2281746
[48,] 43.5218254 1085.6031746
[49,] -288.3531746 43.5218254
[50,] -372.9781746 -288.3531746
[51,] -198.1031746 -372.9781746
[52,] -221.1031746 -198.1031746
[53,] -133.6031746 -221.1031746
[54,] -107.7281746 -133.6031746
[55,] -149.8531746 -107.7281746
[56,] 0.8968254 -149.8531746
[57,] 102.1468254 0.8968254
[58,] -130.9781746 102.1468254
[59,] 144.3968254 -130.9781746
[60,] -221.6845238 144.3968254
[61,] 587.4404762 -221.6845238
[62,] 1190.8154762 587.4404762
[63,] -9.3095238 1190.8154762
[64,] 192.6904762 -9.3095238
[65,] 97.1904762 192.6904762
[66,] -549.9345238 97.1904762
[67,] -424.0595238 -549.9345238
[68,] -173.3095238 -424.0595238
[69,] 37.9404762 -173.3095238
[70,] -225.1845238 37.9404762
[71,] 127.1904762 -225.1845238
[72,] -505.8908730 127.1904762
[73,] -193.7658730 -505.8908730
[74,] 310.6091270 -193.7658730
[75,] -107.5158730 310.6091270
[76,] -7.5158730 -107.5158730
[77,] 104.9841270 -7.5158730
[78,] 984.8591270 104.9841270
[79,] -178.2658730 984.8591270
[80,] 47.4841270 -178.2658730
[81,] -188.2658730 47.4841270
[82,] -231.3908730 -188.2658730
[83,] -487.0158730 -231.3908730
[84,] -675.0972222 -487.0158730
[85,] 391.0277778 -675.0972222
[86,] -87.5972222 391.0277778
[87,] 101.2777778 -87.5972222
[88,] -217.7222222 101.2777778
[89,] -144.2222222 -217.7222222
[90,] -240.3472222 -144.2222222
[91,] -235.4722222 -240.3472222
[92,] 69.2777778 -235.4722222
[93,] 167.5277778 69.2777778
[94,] 624.4027778 167.5277778
[95,] -524.2222222 624.4027778
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 91.4722222 1491.3472222
2 -680.1527778 91.4722222
3 -313.2777778 -680.1527778
4 -278.2777778 -313.2777778
5 -61.7777778 -278.2777778
6 -194.9027778 -61.7777778
7 61.9722222 -194.9027778
8 -204.2777778 61.9722222
9 -141.0277778 -204.2777778
10 82.8472222 -141.0277778
11 -646.7777778 82.8472222
12 -363.8591270 -646.7777778
13 -563.7341270 -363.8591270
14 -303.3591270 -563.7341270
15 -116.4841270 -303.3591270
16 397.5158730 -116.4841270
17 93.0158730 397.5158730
18 -79.1091270 93.0158730
19 137.7658730 -79.1091270
20 -150.4841270 137.7658730
21 -483.2341270 -150.4841270
22 -224.3591270 -483.2341270
23 17.0158730 -224.3591270
24 684.9345238 17.0158730
25 -13.9404762 684.9345238
26 -193.5654762 -13.9404762
27 109.3095238 -193.5654762
28 35.3095238 109.3095238
29 -79.1904762 35.3095238
30 81.6845238 -79.1904762
31 201.5595238 81.6845238
32 -25.6904762 201.5595238
33 210.5595238 -25.6904762
34 -129.5654762 210.5595238
35 283.8095238 -129.5654762
36 -453.2718254 283.8095238
37 -10.1468254 -453.2718254
38 136.2281746 -10.1468254
39 534.1031746 136.2281746
40 99.1031746 534.1031746
41 123.6031746 99.1031746
42 105.4781746 123.6031746
43 586.3531746 105.4781746
44 436.1031746 586.3531746
45 294.3531746 436.1031746
46 234.2281746 294.3531746
47 1085.6031746 234.2281746
48 43.5218254 1085.6031746
49 -288.3531746 43.5218254
50 -372.9781746 -288.3531746
51 -198.1031746 -372.9781746
52 -221.1031746 -198.1031746
53 -133.6031746 -221.1031746
54 -107.7281746 -133.6031746
55 -149.8531746 -107.7281746
56 0.8968254 -149.8531746
57 102.1468254 0.8968254
58 -130.9781746 102.1468254
59 144.3968254 -130.9781746
60 -221.6845238 144.3968254
61 587.4404762 -221.6845238
62 1190.8154762 587.4404762
63 -9.3095238 1190.8154762
64 192.6904762 -9.3095238
65 97.1904762 192.6904762
66 -549.9345238 97.1904762
67 -424.0595238 -549.9345238
68 -173.3095238 -424.0595238
69 37.9404762 -173.3095238
70 -225.1845238 37.9404762
71 127.1904762 -225.1845238
72 -505.8908730 127.1904762
73 -193.7658730 -505.8908730
74 310.6091270 -193.7658730
75 -107.5158730 310.6091270
76 -7.5158730 -107.5158730
77 104.9841270 -7.5158730
78 984.8591270 104.9841270
79 -178.2658730 984.8591270
80 47.4841270 -178.2658730
81 -188.2658730 47.4841270
82 -231.3908730 -188.2658730
83 -487.0158730 -231.3908730
84 -675.0972222 -487.0158730
85 391.0277778 -675.0972222
86 -87.5972222 391.0277778
87 101.2777778 -87.5972222
88 -217.7222222 101.2777778
89 -144.2222222 -217.7222222
90 -240.3472222 -144.2222222
91 -235.4722222 -240.3472222
92 69.2777778 -235.4722222
93 167.5277778 69.2777778
94 624.4027778 167.5277778
95 -524.2222222 624.4027778
> plot(z,main=paste('Residual Lag plot, lowess, and regression line'), ylab='values of Residuals', xlab='lagged values of Residuals')
> lines(lowess(z))
> abline(lm(z))
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/7zc041322577850.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/89ucp1322577850.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/9exk01322577850.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/wessaorg/rcomp/tmp/10hkax1322577850.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/wessaorg/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/wessaorg/rcomp/tmp/11nl2p1322577850.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,hyperlink('http://www.xycoon.com/ols1.htm','Multiple Linear Regression - Ordinary Least Squares',''), 6, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Variable',header=TRUE)
> a<-table.element(a,'Parameter',header=TRUE)
> a<-table.element(a,'S.D.',header=TRUE)
> a<-table.element(a,'T-STAT
H0: parameter = 0',header=TRUE)
> a<-table.element(a,'2-tail p-value',header=TRUE)
> a<-table.element(a,'1-tail p-value',header=TRUE)
> a<-table.row.end(a)
> for (i in 1:k){
+ a<-table.row.start(a)
+ a<-table.element(a,rownames(mysum$coefficients)[i],header=TRUE)
+ a<-table.element(a,mysum$coefficients[i,1])
+ a<-table.element(a, round(mysum$coefficients[i,2],6))
+ a<-table.element(a, round(mysum$coefficients[i,3],4))
+ a<-table.element(a, round(mysum$coefficients[i,4],6))
+ a<-table.element(a, round(mysum$coefficients[i,4]/2,6))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/wessaorg/rcomp/tmp/121ecb1322577850.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Regression Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple R',1,TRUE)
> a<-table.element(a, sqrt(mysum$r.squared))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'R-squared',1,TRUE)
> a<-table.element(a, mysum$r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Adjusted R-squared',1,TRUE)
> a<-table.element(a, mysum$adj.r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (value)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[1])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[2])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[3])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'p-value',1,TRUE)
> a<-table.element(a, 1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3]))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Residual Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Residual Standard Deviation',1,TRUE)
> a<-table.element(a, mysum$sigma)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Sum Squared Residuals',1,TRUE)
> a<-table.element(a, sum(myerror*myerror))
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/wessaorg/rcomp/tmp/13a4bi1322577850.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Actuals, Interpolation, and Residuals', 4, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Time or Index', 1, TRUE)
> a<-table.element(a, 'Actuals', 1, TRUE)
> a<-table.element(a, 'Interpolation
Forecast', 1, TRUE)
> a<-table.element(a, 'Residuals
Prediction Error', 1, TRUE)
> a<-table.row.end(a)
> for (i in 1:n) {
+ a<-table.row.start(a)
+ a<-table.element(a,i, 1, TRUE)
+ a<-table.element(a,x[i])
+ a<-table.element(a,x[i]-mysum$resid[i])
+ a<-table.element(a,mysum$resid[i])
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/wessaorg/rcomp/tmp/14czpc1322577850.tab")
> if (n > n25) {
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'p-values',header=TRUE)
+ a<-table.element(a,'Alternative Hypothesis',3,header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'breakpoint index',header=TRUE)
+ a<-table.element(a,'greater',header=TRUE)
+ a<-table.element(a,'2-sided',header=TRUE)
+ a<-table.element(a,'less',header=TRUE)
+ a<-table.row.end(a)
+ for (mypoint in kp3:nmkm3) {
+ a<-table.row.start(a)
+ a<-table.element(a,mypoint,header=TRUE)
+ a<-table.element(a,gqarr[mypoint-kp3+1,1])
+ a<-table.element(a,gqarr[mypoint-kp3+1,2])
+ a<-table.element(a,gqarr[mypoint-kp3+1,3])
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/var/wessaorg/rcomp/tmp/15bqhl1322577850.tab")
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Meta Analysis of Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'Description',header=TRUE)
+ a<-table.element(a,'# significant tests',header=TRUE)
+ a<-table.element(a,'% significant tests',header=TRUE)
+ a<-table.element(a,'OK/NOK',header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'1% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant1)
+ a<-table.element(a,numsignificant1/numgqtests)
+ if (numsignificant1/numgqtests < 0.01) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'5% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant5)
+ a<-table.element(a,numsignificant5/numgqtests)
+ if (numsignificant5/numgqtests < 0.05) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'10% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant10)
+ a<-table.element(a,numsignificant10/numgqtests)
+ if (numsignificant10/numgqtests < 0.1) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.end(a)
+ table.save(a,file="/var/wessaorg/rcomp/tmp/16ltaz1322577850.tab")
+ }
>
> try(system("convert tmp/1gjku1322577850.ps tmp/1gjku1322577850.png",intern=TRUE))
character(0)
> try(system("convert tmp/2lj0a1322577850.ps tmp/2lj0a1322577850.png",intern=TRUE))
character(0)
> try(system("convert tmp/3f3a51322577850.ps tmp/3f3a51322577850.png",intern=TRUE))
character(0)
> try(system("convert tmp/41gbn1322577850.ps tmp/41gbn1322577850.png",intern=TRUE))
character(0)
> try(system("convert tmp/59wrz1322577850.ps tmp/59wrz1322577850.png",intern=TRUE))
character(0)
> try(system("convert tmp/678n71322577850.ps tmp/678n71322577850.png",intern=TRUE))
character(0)
> try(system("convert tmp/7zc041322577850.ps tmp/7zc041322577850.png",intern=TRUE))
character(0)
> try(system("convert tmp/89ucp1322577850.ps tmp/89ucp1322577850.png",intern=TRUE))
character(0)
> try(system("convert tmp/9exk01322577850.ps tmp/9exk01322577850.png",intern=TRUE))
character(0)
> try(system("convert tmp/10hkax1322577850.ps tmp/10hkax1322577850.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.626 0.496 4.248