R version 2.9.0 (2009-04-17)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
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(332,5140,369,4749,384,3635,373,4305,378,5805,426,4260,423,3869,397,7325,422,9280,409,6222,430,3272,412,7598,470,1345,491,1900,504,1480,484,1472,474,3823,508,4454,492,3357,452,5393,457,8329,457,4152,471,4042,451,7747,493,1451,514,911,522,-406,490,1387,484,2150,506,1577,501,2642,462,4273,465,8064,454,3243,464,1112,427,2280,460,505,473,744,465,-1369,422,-531,415,1041,413,2076,420,577,363,5080,376,6584,380,3761,384,294,346,5020,389,1141,407,3805,393,2127,346,2531,348,3682,353,3263,364,2798,305,5936,307,10568,312,5296,312,1870,286,4390,324,3707,336,5201,327,3748,302,5282,299,5349,311,6249,315,5517,264,8640,278,15767,278,8850,287,5582,279,6496,324,3255,354,6189,354,6452,360,5099,363,6833,385,7046,412,7739,370,10142,389,16054,395,7721,417,6182,404,6490,456,3704,478,6235,468,4655,437,5072,432,3640,441,5147,449,5703,386,11889,396,15603,394,9589),dim=c(2,94),dimnames=list(c('werklooshiedsstotaal','bevolkingstotaal'),1:94))
> y <- array(NA,dim=c(2,94),dimnames=list(c('werklooshiedsstotaal','bevolkingstotaal'),1:94))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No 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
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> 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
werklooshiedsstotaal bevolkingstotaal M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 332 5140 1 0 0 0 0 0 0 0 0 0 0
2 369 4749 0 1 0 0 0 0 0 0 0 0 0
3 384 3635 0 0 1 0 0 0 0 0 0 0 0
4 373 4305 0 0 0 1 0 0 0 0 0 0 0
5 378 5805 0 0 0 0 1 0 0 0 0 0 0
6 426 4260 0 0 0 0 0 1 0 0 0 0 0
7 423 3869 0 0 0 0 0 0 1 0 0 0 0
8 397 7325 0 0 0 0 0 0 0 1 0 0 0
9 422 9280 0 0 0 0 0 0 0 0 1 0 0
10 409 6222 0 0 0 0 0 0 0 0 0 1 0
11 430 3272 0 0 0 0 0 0 0 0 0 0 1
12 412 7598 0 0 0 0 0 0 0 0 0 0 0
13 470 1345 1 0 0 0 0 0 0 0 0 0 0
14 491 1900 0 1 0 0 0 0 0 0 0 0 0
15 504 1480 0 0 1 0 0 0 0 0 0 0 0
16 484 1472 0 0 0 1 0 0 0 0 0 0 0
17 474 3823 0 0 0 0 1 0 0 0 0 0 0
18 508 4454 0 0 0 0 0 1 0 0 0 0 0
19 492 3357 0 0 0 0 0 0 1 0 0 0 0
20 452 5393 0 0 0 0 0 0 0 1 0 0 0
21 457 8329 0 0 0 0 0 0 0 0 1 0 0
22 457 4152 0 0 0 0 0 0 0 0 0 1 0
23 471 4042 0 0 0 0 0 0 0 0 0 0 1
24 451 7747 0 0 0 0 0 0 0 0 0 0 0
25 493 1451 1 0 0 0 0 0 0 0 0 0 0
26 514 911 0 1 0 0 0 0 0 0 0 0 0
27 522 -406 0 0 1 0 0 0 0 0 0 0 0
28 490 1387 0 0 0 1 0 0 0 0 0 0 0
29 484 2150 0 0 0 0 1 0 0 0 0 0 0
30 506 1577 0 0 0 0 0 1 0 0 0 0 0
31 501 2642 0 0 0 0 0 0 1 0 0 0 0
32 462 4273 0 0 0 0 0 0 0 1 0 0 0
33 465 8064 0 0 0 0 0 0 0 0 1 0 0
34 454 3243 0 0 0 0 0 0 0 0 0 1 0
35 464 1112 0 0 0 0 0 0 0 0 0 0 1
36 427 2280 0 0 0 0 0 0 0 0 0 0 0
37 460 505 1 0 0 0 0 0 0 0 0 0 0
38 473 744 0 1 0 0 0 0 0 0 0 0 0
39 465 -1369 0 0 1 0 0 0 0 0 0 0 0
40 422 -531 0 0 0 1 0 0 0 0 0 0 0
41 415 1041 0 0 0 0 1 0 0 0 0 0 0
42 413 2076 0 0 0 0 0 1 0 0 0 0 0
43 420 577 0 0 0 0 0 0 1 0 0 0 0
44 363 5080 0 0 0 0 0 0 0 1 0 0 0
45 376 6584 0 0 0 0 0 0 0 0 1 0 0
46 380 3761 0 0 0 0 0 0 0 0 0 1 0
47 384 294 0 0 0 0 0 0 0 0 0 0 1
48 346 5020 0 0 0 0 0 0 0 0 0 0 0
49 389 1141 1 0 0 0 0 0 0 0 0 0 0
50 407 3805 0 1 0 0 0 0 0 0 0 0 0
51 393 2127 0 0 1 0 0 0 0 0 0 0 0
52 346 2531 0 0 0 1 0 0 0 0 0 0 0
53 348 3682 0 0 0 0 1 0 0 0 0 0 0
54 353 3263 0 0 0 0 0 1 0 0 0 0 0
55 364 2798 0 0 0 0 0 0 1 0 0 0 0
56 305 5936 0 0 0 0 0 0 0 1 0 0 0
57 307 10568 0 0 0 0 0 0 0 0 1 0 0
58 312 5296 0 0 0 0 0 0 0 0 0 1 0
59 312 1870 0 0 0 0 0 0 0 0 0 0 1
60 286 4390 0 0 0 0 0 0 0 0 0 0 0
61 324 3707 1 0 0 0 0 0 0 0 0 0 0
62 336 5201 0 1 0 0 0 0 0 0 0 0 0
63 327 3748 0 0 1 0 0 0 0 0 0 0 0
64 302 5282 0 0 0 1 0 0 0 0 0 0 0
65 299 5349 0 0 0 0 1 0 0 0 0 0 0
66 311 6249 0 0 0 0 0 1 0 0 0 0 0
67 315 5517 0 0 0 0 0 0 1 0 0 0 0
68 264 8640 0 0 0 0 0 0 0 1 0 0 0
69 278 15767 0 0 0 0 0 0 0 0 1 0 0
70 278 8850 0 0 0 0 0 0 0 0 0 1 0
71 287 5582 0 0 0 0 0 0 0 0 0 0 1
72 279 6496 0 0 0 0 0 0 0 0 0 0 0
73 324 3255 1 0 0 0 0 0 0 0 0 0 0
74 354 6189 0 1 0 0 0 0 0 0 0 0 0
75 354 6452 0 0 1 0 0 0 0 0 0 0 0
76 360 5099 0 0 0 1 0 0 0 0 0 0 0
77 363 6833 0 0 0 0 1 0 0 0 0 0 0
78 385 7046 0 0 0 0 0 1 0 0 0 0 0
79 412 7739 0 0 0 0 0 0 1 0 0 0 0
80 370 10142 0 0 0 0 0 0 0 1 0 0 0
81 389 16054 0 0 0 0 0 0 0 0 1 0 0
82 395 7721 0 0 0 0 0 0 0 0 0 1 0
83 417 6182 0 0 0 0 0 0 0 0 0 0 1
84 404 6490 0 0 0 0 0 0 0 0 0 0 0
85 456 3704 1 0 0 0 0 0 0 0 0 0 0
86 478 6235 0 1 0 0 0 0 0 0 0 0 0
87 468 4655 0 0 1 0 0 0 0 0 0 0 0
88 437 5072 0 0 0 1 0 0 0 0 0 0 0
89 432 3640 0 0 0 0 1 0 0 0 0 0 0
90 441 5147 0 0 0 0 0 1 0 0 0 0 0
91 449 5703 0 0 0 0 0 0 1 0 0 0 0
92 386 11889 0 0 0 0 0 0 0 1 0 0 0
93 396 15603 0 0 0 0 0 0 0 0 1 0 0
94 394 9589 0 0 0 0 0 0 0 0 0 1 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) bevolkingstotaal M1 M2
440.69616 -0.01199 -4.34812 31.61965
M3 M4 M5 M6
16.88779 -2.04979 6.87509 28.24652
M7 M8 M9 M10
29.56873 22.12640 80.82054 17.37205
M11
-7.40525
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-102.058 -46.221 8.871 52.183 103.194
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 440.69616 28.40997 15.512 < 2e-16 ***
bevolkingstotaal -0.01199 0.00284 -4.223 6.27e-05 ***
M1 -4.34812 33.18225 -0.131 0.8961
M2 31.61965 32.42643 0.975 0.3324
M3 16.88779 33.17509 0.509 0.6121
M4 -2.04979 32.79338 -0.063 0.9503
M5 6.87509 32.27813 0.213 0.8319
M6 28.24652 32.19241 0.877 0.3828
M7 29.56873 32.28449 0.916 0.3624
M8 22.12640 32.25366 0.686 0.4947
M9 80.82054 35.62020 2.269 0.0259 *
M10 17.37205 31.94389 0.544 0.5881
M11 -7.40525 33.74182 -0.219 0.8268
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 61.68 on 81 degrees of freedom
Multiple R-squared: 0.2457, Adjusted R-squared: 0.134
F-statistic: 2.199 on 12 and 81 DF, p-value: 0.01911
> 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.0159462601 0.0318925202 0.984053740
[2,] 0.0039379862 0.0078759723 0.996062014
[3,] 0.0514678023 0.1029356045 0.948532198
[4,] 0.0369035106 0.0738070212 0.963096489
[5,] 0.0187548728 0.0375097456 0.981245127
[6,] 0.0078279953 0.0156559906 0.992172005
[7,] 0.0044790225 0.0089580451 0.995520977
[8,] 0.0072146990 0.0144293979 0.992785301
[9,] 0.0063082822 0.0126165644 0.993691718
[10,] 0.0037833156 0.0075666312 0.996216684
[11,] 0.0018525554 0.0037051109 0.998147445
[12,] 0.0013116423 0.0026232847 0.998688358
[13,] 0.0007059416 0.0014118832 0.999294058
[14,] 0.0006079808 0.0012159616 0.999392019
[15,] 0.0009051064 0.0018102128 0.999094894
[16,] 0.0005570134 0.0011140267 0.999442987
[17,] 0.0004466583 0.0008933166 0.999553342
[18,] 0.0002592700 0.0005185399 0.999740730
[19,] 0.0002534240 0.0005068481 0.999746576
[20,] 0.0005707006 0.0011414012 0.999429299
[21,] 0.0153990380 0.0307980760 0.984600962
[22,] 0.0121237561 0.0242475122 0.987876244
[23,] 0.0092662702 0.0185325404 0.990733730
[24,] 0.0122850374 0.0245700748 0.987714963
[25,] 0.0198506136 0.0397012272 0.980149386
[26,] 0.0279403436 0.0558806872 0.972059656
[27,] 0.0408544541 0.0817089081 0.959145546
[28,] 0.0565219623 0.1130439246 0.943478038
[29,] 0.0663262317 0.1326524634 0.933673768
[30,] 0.0925522336 0.1851044672 0.907447766
[31,] 0.0963479346 0.1926958691 0.903652065
[32,] 0.1319007982 0.2638015963 0.868099202
[33,] 0.1455002309 0.2910004617 0.854499769
[34,] 0.1334205321 0.2668410643 0.866579468
[35,] 0.1109223496 0.2218446991 0.889077650
[36,] 0.1019656354 0.2039312707 0.898034365
[37,] 0.1037668345 0.2075336690 0.896233166
[38,] 0.1016258911 0.2032517823 0.898374109
[39,] 0.1150282846 0.2300565693 0.884971715
[40,] 0.1137492376 0.2274984753 0.886250762
[41,] 0.1215913164 0.2431826329 0.878408684
[42,] 0.1158201381 0.2316402762 0.884179862
[43,] 0.1129402621 0.2258805242 0.887059738
[44,] 0.1261374226 0.2522748452 0.873862577
[45,] 0.1476323030 0.2952646060 0.852367697
[46,] 0.1429848480 0.2859696960 0.857015152
[47,] 0.1382765436 0.2765530872 0.861723456
[48,] 0.1401797839 0.2803595678 0.859820216
[49,] 0.1457104659 0.2914209319 0.854289534
[50,] 0.1617958773 0.3235917547 0.838204123
[51,] 0.1780230626 0.3560461253 0.821976937
[52,] 0.2140759068 0.4281518136 0.785924093
[53,] 0.3028104428 0.6056208857 0.697189557
[54,] 0.3551119606 0.7102239212 0.644888039
[55,] 0.4221546113 0.8443092225 0.577845389
[56,] 0.5442368705 0.9115262590 0.455763129
[57,] 0.6560356589 0.6879286823 0.343964341
[58,] 0.8331782005 0.3336435990 0.166821800
[59,] 0.9465598195 0.1068803610 0.053440181
[60,] 0.9696758557 0.0606482885 0.030324144
[61,] 0.9910596778 0.0178806443 0.008940322
[62,] 0.9795535922 0.0408928155 0.020446408
[63,] 0.9772425416 0.0455149168 0.022757458
> postscript(file="/var/www/html/rcomp/tmp/1d4v81293373499.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/www/html/rcomp/tmp/2d4v81293373499.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/www/html/rcomp/tmp/3ovct1293373499.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/www/html/rcomp/tmp/4ovct1293373499.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/www/html/rcomp/tmp/5ovct1293373499.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 = 94
Frequency = 1
1 2 3 4 5
-42.71670026 -46.37277265 -29.99836436 -14.02712469 0.03379372
6 7 8 9 10
8.13699053 -0.87352073 22.00809231 11.75544531 25.53688522
11 12 13 14 15
35.94211328 62.40790895 49.77922825 41.46619977 64.16203773
16 17 18 19 20
63.00369624 72.26855796 92.46315387 61.98732654 53.84238318
21 22 23 24 25
35.35244873 48.71648259 86.17482344 103.19449832 74.05022471
26 27 28 29 30
52.60756295 59.54789311 67.98450097 62.20839680 55.96639137
31 32 33 34 35
62.41409568 50.41298659 40.17495758 34.81708839 44.04256271
36 37 38 39 40
13.64225619 29.70718081 9.60514400 -8.99898985 -23.01334070
41 42 43 44 45
-20.08910394 -31.05033282 -43.34635429 -38.91065354 -66.57103078
46 47 48 49 50
-32.97181569 -45.76569302 -34.50368428 -33.66684042 -19.69183549
51 52 53 54 55
-39.08008763 -62.29832966 -55.42210715 -76.81757054 -72.71538122
56 57 58 59 60
-86.64675758 -87.80074861 -82.56634803 -98.86861353 -102.05771987
61 62 63 64 65
-67.89913358 -73.95305188 -85.64343416 -73.31237427 -84.43388918
66 67 68 69 70
-83.01383998 -89.11312289 -95.22435723 -54.46196925 -73.95199491
71 72 73 74 75
-79.35975624 -83.80565806 -73.31885434 -44.10640560 -26.22103382
76 77 78 79 80
-17.50664175 -2.63993869 0.54261456 34.52984071 28.78542285
81 82 83 84 85
59.97931363 29.51069370 57.83456336 41.12239875 64.06489483
86 87 88 89 90
80.44515890 66.23197897 59.16961387 28.07429048 33.77259301
91 92 93 94
47.11711619 65.73288342 61.57158339 50.90900873
> postscript(file="/var/www/html/rcomp/tmp/6gncw1293373499.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 = 94
Frequency = 1
lag(myerror, k = 1) myerror
0 -42.71670026 NA
1 -46.37277265 -42.71670026
2 -29.99836436 -46.37277265
3 -14.02712469 -29.99836436
4 0.03379372 -14.02712469
5 8.13699053 0.03379372
6 -0.87352073 8.13699053
7 22.00809231 -0.87352073
8 11.75544531 22.00809231
9 25.53688522 11.75544531
10 35.94211328 25.53688522
11 62.40790895 35.94211328
12 49.77922825 62.40790895
13 41.46619977 49.77922825
14 64.16203773 41.46619977
15 63.00369624 64.16203773
16 72.26855796 63.00369624
17 92.46315387 72.26855796
18 61.98732654 92.46315387
19 53.84238318 61.98732654
20 35.35244873 53.84238318
21 48.71648259 35.35244873
22 86.17482344 48.71648259
23 103.19449832 86.17482344
24 74.05022471 103.19449832
25 52.60756295 74.05022471
26 59.54789311 52.60756295
27 67.98450097 59.54789311
28 62.20839680 67.98450097
29 55.96639137 62.20839680
30 62.41409568 55.96639137
31 50.41298659 62.41409568
32 40.17495758 50.41298659
33 34.81708839 40.17495758
34 44.04256271 34.81708839
35 13.64225619 44.04256271
36 29.70718081 13.64225619
37 9.60514400 29.70718081
38 -8.99898985 9.60514400
39 -23.01334070 -8.99898985
40 -20.08910394 -23.01334070
41 -31.05033282 -20.08910394
42 -43.34635429 -31.05033282
43 -38.91065354 -43.34635429
44 -66.57103078 -38.91065354
45 -32.97181569 -66.57103078
46 -45.76569302 -32.97181569
47 -34.50368428 -45.76569302
48 -33.66684042 -34.50368428
49 -19.69183549 -33.66684042
50 -39.08008763 -19.69183549
51 -62.29832966 -39.08008763
52 -55.42210715 -62.29832966
53 -76.81757054 -55.42210715
54 -72.71538122 -76.81757054
55 -86.64675758 -72.71538122
56 -87.80074861 -86.64675758
57 -82.56634803 -87.80074861
58 -98.86861353 -82.56634803
59 -102.05771987 -98.86861353
60 -67.89913358 -102.05771987
61 -73.95305188 -67.89913358
62 -85.64343416 -73.95305188
63 -73.31237427 -85.64343416
64 -84.43388918 -73.31237427
65 -83.01383998 -84.43388918
66 -89.11312289 -83.01383998
67 -95.22435723 -89.11312289
68 -54.46196925 -95.22435723
69 -73.95199491 -54.46196925
70 -79.35975624 -73.95199491
71 -83.80565806 -79.35975624
72 -73.31885434 -83.80565806
73 -44.10640560 -73.31885434
74 -26.22103382 -44.10640560
75 -17.50664175 -26.22103382
76 -2.63993869 -17.50664175
77 0.54261456 -2.63993869
78 34.52984071 0.54261456
79 28.78542285 34.52984071
80 59.97931363 28.78542285
81 29.51069370 59.97931363
82 57.83456336 29.51069370
83 41.12239875 57.83456336
84 64.06489483 41.12239875
85 80.44515890 64.06489483
86 66.23197897 80.44515890
87 59.16961387 66.23197897
88 28.07429048 59.16961387
89 33.77259301 28.07429048
90 47.11711619 33.77259301
91 65.73288342 47.11711619
92 61.57158339 65.73288342
93 50.90900873 61.57158339
94 NA 50.90900873
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -46.37277265 -42.71670026
[2,] -29.99836436 -46.37277265
[3,] -14.02712469 -29.99836436
[4,] 0.03379372 -14.02712469
[5,] 8.13699053 0.03379372
[6,] -0.87352073 8.13699053
[7,] 22.00809231 -0.87352073
[8,] 11.75544531 22.00809231
[9,] 25.53688522 11.75544531
[10,] 35.94211328 25.53688522
[11,] 62.40790895 35.94211328
[12,] 49.77922825 62.40790895
[13,] 41.46619977 49.77922825
[14,] 64.16203773 41.46619977
[15,] 63.00369624 64.16203773
[16,] 72.26855796 63.00369624
[17,] 92.46315387 72.26855796
[18,] 61.98732654 92.46315387
[19,] 53.84238318 61.98732654
[20,] 35.35244873 53.84238318
[21,] 48.71648259 35.35244873
[22,] 86.17482344 48.71648259
[23,] 103.19449832 86.17482344
[24,] 74.05022471 103.19449832
[25,] 52.60756295 74.05022471
[26,] 59.54789311 52.60756295
[27,] 67.98450097 59.54789311
[28,] 62.20839680 67.98450097
[29,] 55.96639137 62.20839680
[30,] 62.41409568 55.96639137
[31,] 50.41298659 62.41409568
[32,] 40.17495758 50.41298659
[33,] 34.81708839 40.17495758
[34,] 44.04256271 34.81708839
[35,] 13.64225619 44.04256271
[36,] 29.70718081 13.64225619
[37,] 9.60514400 29.70718081
[38,] -8.99898985 9.60514400
[39,] -23.01334070 -8.99898985
[40,] -20.08910394 -23.01334070
[41,] -31.05033282 -20.08910394
[42,] -43.34635429 -31.05033282
[43,] -38.91065354 -43.34635429
[44,] -66.57103078 -38.91065354
[45,] -32.97181569 -66.57103078
[46,] -45.76569302 -32.97181569
[47,] -34.50368428 -45.76569302
[48,] -33.66684042 -34.50368428
[49,] -19.69183549 -33.66684042
[50,] -39.08008763 -19.69183549
[51,] -62.29832966 -39.08008763
[52,] -55.42210715 -62.29832966
[53,] -76.81757054 -55.42210715
[54,] -72.71538122 -76.81757054
[55,] -86.64675758 -72.71538122
[56,] -87.80074861 -86.64675758
[57,] -82.56634803 -87.80074861
[58,] -98.86861353 -82.56634803
[59,] -102.05771987 -98.86861353
[60,] -67.89913358 -102.05771987
[61,] -73.95305188 -67.89913358
[62,] -85.64343416 -73.95305188
[63,] -73.31237427 -85.64343416
[64,] -84.43388918 -73.31237427
[65,] -83.01383998 -84.43388918
[66,] -89.11312289 -83.01383998
[67,] -95.22435723 -89.11312289
[68,] -54.46196925 -95.22435723
[69,] -73.95199491 -54.46196925
[70,] -79.35975624 -73.95199491
[71,] -83.80565806 -79.35975624
[72,] -73.31885434 -83.80565806
[73,] -44.10640560 -73.31885434
[74,] -26.22103382 -44.10640560
[75,] -17.50664175 -26.22103382
[76,] -2.63993869 -17.50664175
[77,] 0.54261456 -2.63993869
[78,] 34.52984071 0.54261456
[79,] 28.78542285 34.52984071
[80,] 59.97931363 28.78542285
[81,] 29.51069370 59.97931363
[82,] 57.83456336 29.51069370
[83,] 41.12239875 57.83456336
[84,] 64.06489483 41.12239875
[85,] 80.44515890 64.06489483
[86,] 66.23197897 80.44515890
[87,] 59.16961387 66.23197897
[88,] 28.07429048 59.16961387
[89,] 33.77259301 28.07429048
[90,] 47.11711619 33.77259301
[91,] 65.73288342 47.11711619
[92,] 61.57158339 65.73288342
[93,] 50.90900873 61.57158339
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -46.37277265 -42.71670026
2 -29.99836436 -46.37277265
3 -14.02712469 -29.99836436
4 0.03379372 -14.02712469
5 8.13699053 0.03379372
6 -0.87352073 8.13699053
7 22.00809231 -0.87352073
8 11.75544531 22.00809231
9 25.53688522 11.75544531
10 35.94211328 25.53688522
11 62.40790895 35.94211328
12 49.77922825 62.40790895
13 41.46619977 49.77922825
14 64.16203773 41.46619977
15 63.00369624 64.16203773
16 72.26855796 63.00369624
17 92.46315387 72.26855796
18 61.98732654 92.46315387
19 53.84238318 61.98732654
20 35.35244873 53.84238318
21 48.71648259 35.35244873
22 86.17482344 48.71648259
23 103.19449832 86.17482344
24 74.05022471 103.19449832
25 52.60756295 74.05022471
26 59.54789311 52.60756295
27 67.98450097 59.54789311
28 62.20839680 67.98450097
29 55.96639137 62.20839680
30 62.41409568 55.96639137
31 50.41298659 62.41409568
32 40.17495758 50.41298659
33 34.81708839 40.17495758
34 44.04256271 34.81708839
35 13.64225619 44.04256271
36 29.70718081 13.64225619
37 9.60514400 29.70718081
38 -8.99898985 9.60514400
39 -23.01334070 -8.99898985
40 -20.08910394 -23.01334070
41 -31.05033282 -20.08910394
42 -43.34635429 -31.05033282
43 -38.91065354 -43.34635429
44 -66.57103078 -38.91065354
45 -32.97181569 -66.57103078
46 -45.76569302 -32.97181569
47 -34.50368428 -45.76569302
48 -33.66684042 -34.50368428
49 -19.69183549 -33.66684042
50 -39.08008763 -19.69183549
51 -62.29832966 -39.08008763
52 -55.42210715 -62.29832966
53 -76.81757054 -55.42210715
54 -72.71538122 -76.81757054
55 -86.64675758 -72.71538122
56 -87.80074861 -86.64675758
57 -82.56634803 -87.80074861
58 -98.86861353 -82.56634803
59 -102.05771987 -98.86861353
60 -67.89913358 -102.05771987
61 -73.95305188 -67.89913358
62 -85.64343416 -73.95305188
63 -73.31237427 -85.64343416
64 -84.43388918 -73.31237427
65 -83.01383998 -84.43388918
66 -89.11312289 -83.01383998
67 -95.22435723 -89.11312289
68 -54.46196925 -95.22435723
69 -73.95199491 -54.46196925
70 -79.35975624 -73.95199491
71 -83.80565806 -79.35975624
72 -73.31885434 -83.80565806
73 -44.10640560 -73.31885434
74 -26.22103382 -44.10640560
75 -17.50664175 -26.22103382
76 -2.63993869 -17.50664175
77 0.54261456 -2.63993869
78 34.52984071 0.54261456
79 28.78542285 34.52984071
80 59.97931363 28.78542285
81 29.51069370 59.97931363
82 57.83456336 29.51069370
83 41.12239875 57.83456336
84 64.06489483 41.12239875
85 80.44515890 64.06489483
86 66.23197897 80.44515890
87 59.16961387 66.23197897
88 28.07429048 59.16961387
89 33.77259301 28.07429048
90 47.11711619 33.77259301
91 65.73288342 47.11711619
92 61.57158339 65.73288342
93 50.90900873 61.57158339
> 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/www/html/rcomp/tmp/7rwth1293373499.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/www/html/rcomp/tmp/8rwth1293373499.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/www/html/rcomp/tmp/9rwth1293373499.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/www/html/rcomp/tmp/102nsk1293373499.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/www/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/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/www/html/rcomp/tmp/11no881293373499.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/www/html/rcomp/tmp/129opw1293373499.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/www/html/rcomp/tmp/135g5n1293373499.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/www/html/rcomp/tmp/14qzla1293373499.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/www/html/rcomp/tmp/15tz2g1293373499.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/www/html/rcomp/tmp/16ps3h1293373500.tab")
+ }
>
> try(system("convert tmp/1d4v81293373499.ps tmp/1d4v81293373499.png",intern=TRUE))
character(0)
> try(system("convert tmp/2d4v81293373499.ps tmp/2d4v81293373499.png",intern=TRUE))
character(0)
> try(system("convert tmp/3ovct1293373499.ps tmp/3ovct1293373499.png",intern=TRUE))
character(0)
> try(system("convert tmp/4ovct1293373499.ps tmp/4ovct1293373499.png",intern=TRUE))
character(0)
> try(system("convert tmp/5ovct1293373499.ps tmp/5ovct1293373499.png",intern=TRUE))
character(0)
> try(system("convert tmp/6gncw1293373499.ps tmp/6gncw1293373499.png",intern=TRUE))
character(0)
> try(system("convert tmp/7rwth1293373499.ps tmp/7rwth1293373499.png",intern=TRUE))
character(0)
> try(system("convert tmp/8rwth1293373499.ps tmp/8rwth1293373499.png",intern=TRUE))
character(0)
> try(system("convert tmp/9rwth1293373499.ps tmp/9rwth1293373499.png",intern=TRUE))
character(0)
> try(system("convert tmp/102nsk1293373499.ps tmp/102nsk1293373499.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.962 1.714 7.476