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(370.861
+ ,378.205
+ ,377.632
+ ,369.167
+ ,370.861
+ ,378.205
+ ,371.551
+ ,369.167
+ ,370.861
+ ,382.842
+ ,371.551
+ ,369.167
+ ,381.903
+ ,382.842
+ ,371.551
+ ,384.502
+ ,381.903
+ ,382.842
+ ,392.058
+ ,384.502
+ ,381.903
+ ,384.359
+ ,392.058
+ ,384.502
+ ,388.884
+ ,384.359
+ ,392.058
+ ,386.586
+ ,388.884
+ ,384.359
+ ,387.495
+ ,386.586
+ ,388.884
+ ,385.705
+ ,387.495
+ ,386.586
+ ,378.67
+ ,385.705
+ ,387.495
+ ,377.367
+ ,378.67
+ ,385.705
+ ,376.911
+ ,377.367
+ ,378.67
+ ,389.827
+ ,376.911
+ ,377.367
+ ,387.82
+ ,389.827
+ ,376.911
+ ,387.267
+ ,387.82
+ ,389.827
+ ,380.575
+ ,387.267
+ ,387.82
+ ,372.402
+ ,380.575
+ ,387.267
+ ,376.74
+ ,372.402
+ ,380.575
+ ,377.795
+ ,376.74
+ ,372.402
+ ,376.126
+ ,377.795
+ ,376.74
+ ,370.804
+ ,376.126
+ ,377.795
+ ,367.98
+ ,370.804
+ ,376.126
+ ,367.866
+ ,367.98
+ ,370.804
+ ,366.121
+ ,367.866
+ ,367.98
+ ,379.421
+ ,366.121
+ ,367.866
+ ,378.519
+ ,379.421
+ ,366.121
+ ,372.423
+ ,378.519
+ ,379.421
+ ,355.072
+ ,372.423
+ ,378.519
+ ,344.693
+ ,355.072
+ ,372.423
+ ,342.892
+ ,344.693
+ ,355.072
+ ,344.178
+ ,342.892
+ ,344.693
+ ,337.606
+ ,344.178
+ ,342.892
+ ,327.103
+ ,337.606
+ ,344.178
+ ,323.953
+ ,327.103
+ ,337.606
+ ,316.532
+ ,323.953
+ ,327.103
+ ,306.307
+ ,316.532
+ ,323.953
+ ,327.225
+ ,306.307
+ ,316.532
+ ,329.573
+ ,327.225
+ ,306.307
+ ,313.761
+ ,329.573
+ ,327.225
+ ,307.836
+ ,313.761
+ ,329.573
+ ,300.074
+ ,307.836
+ ,313.761
+ ,304.198
+ ,300.074
+ ,307.836
+ ,306.122
+ ,304.198
+ ,300.074
+ ,300.414
+ ,306.122
+ ,304.198
+ ,292.133
+ ,300.414
+ ,306.122
+ ,290.616
+ ,292.133
+ ,300.414
+ ,280.244
+ ,290.616
+ ,292.133
+ ,285.179
+ ,280.244
+ ,290.616
+ ,305.486
+ ,285.179
+ ,280.244
+ ,305.957
+ ,305.486
+ ,285.179
+ ,293.886
+ ,305.957
+ ,305.486
+ ,289.441
+ ,293.886
+ ,305.957
+ ,288.776
+ ,289.441
+ ,293.886
+ ,299.149
+ ,288.776
+ ,289.441
+ ,306.532
+ ,299.149
+ ,288.776
+ ,309.914
+ ,306.532
+ ,299.149
+ ,313.468
+ ,309.914
+ ,306.532
+ ,314.901
+ ,313.468
+ ,309.914
+ ,309.16
+ ,314.901
+ ,313.468
+ ,316.15
+ ,309.16
+ ,314.901
+ ,336.544
+ ,316.15
+ ,309.16
+ ,339.196
+ ,336.544
+ ,316.15
+ ,326.738
+ ,339.196
+ ,336.544
+ ,320.838
+ ,326.738
+ ,339.196
+ ,318.62
+ ,320.838
+ ,326.738
+ ,331.533
+ ,318.62
+ ,320.838
+ ,335.378
+ ,331.533
+ ,318.62)
+ ,dim=c(3
+ ,70)
+ ,dimnames=list(c('Maandelijksewerkloosheid'
+ ,'y-1'
+ ,'y-2')
+ ,1:70))
> y <- array(NA,dim=c(3,70),dimnames=list(c('Maandelijksewerkloosheid','y-1','y-2'),1:70))
> 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
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
Maandelijksewerkloosheid y-1 y-2 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 370.861 378.205 377.632 1 0 0 0 0 0 0 0 0 0 0
2 369.167 370.861 378.205 0 1 0 0 0 0 0 0 0 0 0
3 371.551 369.167 370.861 0 0 1 0 0 0 0 0 0 0 0
4 382.842 371.551 369.167 0 0 0 1 0 0 0 0 0 0 0
5 381.903 382.842 371.551 0 0 0 0 1 0 0 0 0 0 0
6 384.502 381.903 382.842 0 0 0 0 0 1 0 0 0 0 0
7 392.058 384.502 381.903 0 0 0 0 0 0 1 0 0 0 0
8 384.359 392.058 384.502 0 0 0 0 0 0 0 1 0 0 0
9 388.884 384.359 392.058 0 0 0 0 0 0 0 0 1 0 0
10 386.586 388.884 384.359 0 0 0 0 0 0 0 0 0 1 0
11 387.495 386.586 388.884 0 0 0 0 0 0 0 0 0 0 1
12 385.705 387.495 386.586 0 0 0 0 0 0 0 0 0 0 0
13 378.670 385.705 387.495 1 0 0 0 0 0 0 0 0 0 0
14 377.367 378.670 385.705 0 1 0 0 0 0 0 0 0 0 0
15 376.911 377.367 378.670 0 0 1 0 0 0 0 0 0 0 0
16 389.827 376.911 377.367 0 0 0 1 0 0 0 0 0 0 0
17 387.820 389.827 376.911 0 0 0 0 1 0 0 0 0 0 0
18 387.267 387.820 389.827 0 0 0 0 0 1 0 0 0 0 0
19 380.575 387.267 387.820 0 0 0 0 0 0 1 0 0 0 0
20 372.402 380.575 387.267 0 0 0 0 0 0 0 1 0 0 0
21 376.740 372.402 380.575 0 0 0 0 0 0 0 0 1 0 0
22 377.795 376.740 372.402 0 0 0 0 0 0 0 0 0 1 0
23 376.126 377.795 376.740 0 0 0 0 0 0 0 0 0 0 1
24 370.804 376.126 377.795 0 0 0 0 0 0 0 0 0 0 0
25 367.980 370.804 376.126 1 0 0 0 0 0 0 0 0 0 0
26 367.866 367.980 370.804 0 1 0 0 0 0 0 0 0 0 0
27 366.121 367.866 367.980 0 0 1 0 0 0 0 0 0 0 0
28 379.421 366.121 367.866 0 0 0 1 0 0 0 0 0 0 0
29 378.519 379.421 366.121 0 0 0 0 1 0 0 0 0 0 0
30 372.423 378.519 379.421 0 0 0 0 0 1 0 0 0 0 0
31 355.072 372.423 378.519 0 0 0 0 0 0 1 0 0 0 0
32 344.693 355.072 372.423 0 0 0 0 0 0 0 1 0 0 0
33 342.892 344.693 355.072 0 0 0 0 0 0 0 0 1 0 0
34 344.178 342.892 344.693 0 0 0 0 0 0 0 0 0 1 0
35 337.606 344.178 342.892 0 0 0 0 0 0 0 0 0 0 1
36 327.103 337.606 344.178 0 0 0 0 0 0 0 0 0 0 0
37 323.953 327.103 337.606 1 0 0 0 0 0 0 0 0 0 0
38 316.532 323.953 327.103 0 1 0 0 0 0 0 0 0 0 0
39 306.307 316.532 323.953 0 0 1 0 0 0 0 0 0 0 0
40 327.225 306.307 316.532 0 0 0 1 0 0 0 0 0 0 0
41 329.573 327.225 306.307 0 0 0 0 1 0 0 0 0 0 0
42 313.761 329.573 327.225 0 0 0 0 0 1 0 0 0 0 0
43 307.836 313.761 329.573 0 0 0 0 0 0 1 0 0 0 0
44 300.074 307.836 313.761 0 0 0 0 0 0 0 1 0 0 0
45 304.198 300.074 307.836 0 0 0 0 0 0 0 0 1 0 0
46 306.122 304.198 300.074 0 0 0 0 0 0 0 0 0 1 0
47 300.414 306.122 304.198 0 0 0 0 0 0 0 0 0 0 1
48 292.133 300.414 306.122 0 0 0 0 0 0 0 0 0 0 0
49 290.616 292.133 300.414 1 0 0 0 0 0 0 0 0 0 0
50 280.244 290.616 292.133 0 1 0 0 0 0 0 0 0 0 0
51 285.179 280.244 290.616 0 0 1 0 0 0 0 0 0 0 0
52 305.486 285.179 280.244 0 0 0 1 0 0 0 0 0 0 0
53 305.957 305.486 285.179 0 0 0 0 1 0 0 0 0 0 0
54 293.886 305.957 305.486 0 0 0 0 0 1 0 0 0 0 0
55 289.441 293.886 305.957 0 0 0 0 0 0 1 0 0 0 0
56 288.776 289.441 293.886 0 0 0 0 0 0 0 1 0 0 0
57 299.149 288.776 289.441 0 0 0 0 0 0 0 0 1 0 0
58 306.532 299.149 288.776 0 0 0 0 0 0 0 0 0 1 0
59 309.914 306.532 299.149 0 0 0 0 0 0 0 0 0 0 1
60 313.468 309.914 306.532 0 0 0 0 0 0 0 0 0 0 0
61 314.901 313.468 309.914 1 0 0 0 0 0 0 0 0 0 0
62 309.160 314.901 313.468 0 1 0 0 0 0 0 0 0 0 0
63 316.150 309.160 314.901 0 0 1 0 0 0 0 0 0 0 0
64 336.544 316.150 309.160 0 0 0 1 0 0 0 0 0 0 0
65 339.196 336.544 316.150 0 0 0 0 1 0 0 0 0 0 0
66 326.738 339.196 336.544 0 0 0 0 0 1 0 0 0 0 0
67 320.838 326.738 339.196 0 0 0 0 0 0 1 0 0 0 0
68 318.620 320.838 326.738 0 0 0 0 0 0 0 1 0 0 0
69 331.533 318.620 320.838 0 0 0 0 0 0 0 0 1 0 0
70 335.378 331.533 318.620 0 0 0 0 0 0 0 0 0 1 0
t
1 1
2 2
3 3
4 4
5 5
6 6
7 7
8 8
9 9
10 10
11 11
12 12
13 13
14 14
15 15
16 16
17 17
18 18
19 19
20 20
21 21
22 22
23 23
24 24
25 25
26 26
27 27
28 28
29 29
30 30
31 31
32 32
33 33
34 34
35 35
36 36
37 37
38 38
39 39
40 40
41 41
42 42
43 43
44 44
45 45
46 46
47 47
48 48
49 49
50 50
51 51
52 52
53 53
54 54
55 55
56 56
57 57
58 58
59 59
60 60
61 61
62 62
63 63
64 64
65 65
66 66
67 67
68 68
69 69
70 70
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) `y-1` `y-2` M1 M2 M3
-0.18365 1.20254 -0.21487 1.50327 0.36918 5.28169
M4 M5 M6 M7 M8 M9
20.46136 0.92236 -3.26128 0.22462 -0.95888 10.99847
M10 M11 t
4.95759 1.75349 0.00973
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-11.79254 -2.90994 -0.04397 2.68501 11.62859
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.18365 13.97031 -0.013 0.989559
`y-1` 1.20254 0.13046 9.218 9.51e-13 ***
`y-2` -0.21487 0.13429 -1.600 0.115321
M1 1.50327 3.04816 0.493 0.623854
M2 0.36918 3.04914 0.121 0.904071
M3 5.28169 3.06457 1.723 0.090422 .
M4 20.46136 3.07118 6.662 1.35e-08 ***
M5 0.92236 3.88402 0.237 0.813169
M6 -3.26128 3.06597 -1.064 0.292110
M7 0.22462 3.14022 0.072 0.943236
M8 -0.95888 3.07596 -0.312 0.756420
M9 10.99847 3.08770 3.562 0.000769 ***
M10 4.95759 3.19533 1.552 0.126515
M11 1.75349 3.21088 0.546 0.587198
t 0.00973 0.06012 0.162 0.872023
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.012 on 55 degrees of freedom
Multiple R-squared: 0.9838, Adjusted R-squared: 0.9797
F-statistic: 239.1 on 14 and 55 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.0613928 0.1227855966 0.9386072017
[2,] 0.2375950 0.4751900927 0.7624049537
[3,] 0.5027684 0.9944631405 0.4972315702
[4,] 0.4317350 0.8634699681 0.5682650159
[5,] 0.3406967 0.6813933054 0.6593033473
[6,] 0.2739408 0.5478816507 0.7260591747
[7,] 0.2086027 0.4172054324 0.7913972838
[8,] 0.2189605 0.4379209993 0.7810395004
[9,] 0.2697256 0.5394512312 0.7302743844
[10,] 0.2187306 0.4374612202 0.7812693899
[11,] 0.1583921 0.3167841410 0.8416079295
[12,] 0.1134976 0.2269951277 0.8865024361
[13,] 0.4690175 0.9380349040 0.5309825480
[14,] 0.8819876 0.2360247714 0.1180123857
[15,] 0.8442486 0.3115027904 0.1557513952
[16,] 0.8234735 0.3530529179 0.1765264589
[17,] 0.8111209 0.3777582298 0.1888791149
[18,] 0.7745015 0.4509969749 0.2254984874
[19,] 0.7222548 0.5554904495 0.2777452248
[20,] 0.7238309 0.5523381419 0.2761690709
[21,] 0.7659027 0.4681945828 0.2340972914
[22,] 0.9304273 0.1391454026 0.0695727013
[23,] 0.9938543 0.0122914401 0.0061457201
[24,] 0.9945678 0.0108643341 0.0054321670
[25,] 0.9950568 0.0098863531 0.0049431765
[26,] 0.9980367 0.0039266413 0.0019633206
[27,] 0.9957123 0.0085753612 0.0042876806
[28,] 0.9912510 0.0174980479 0.0087490240
[29,] 0.9991426 0.0017148382 0.0008574191
[30,] 0.9998879 0.0002242423 0.0001121211
[31,] 0.9995348 0.0009304496 0.0004652248
[32,] 0.9987983 0.0024033965 0.0012016983
[33,] 0.9963876 0.0072247761 0.0036123880
[34,] 0.9873158 0.0253683362 0.0126841681
[35,] 0.9559750 0.0880499881 0.0440249941
> postscript(file="/var/www/html/rcomp/tmp/1sli51290976213.ps",horizontal=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/2sli51290976213.ps",horizontal=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/32c0q1290976213.ps",horizontal=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/42c0q1290976213.ps",horizontal=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/52c0q1290976213.ps",horizontal=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 = 70
Frequency = 1
1 2 3 4 5 6
-4.13398365 4.25096541 2.17182653 -4.95742384 0.56718862 10.89539276
7 8 9 10 11 12
11.62859205 -3.42459841 0.01524952 -3.34739188 4.49170561 2.85858576
13 14 15 16 17 18
-3.34154645 4.55507904 -0.76785897 -2.77287537 -0.88062501 7.92905807
19 20 21 22 23 24
-2.02481099 -1.09545118 -0.33406741 -0.22068185 0.96811105 -0.37639566
25 26 27 28 29 30
1.32791596 4.59071519 -2.54622747 -2.36168647 -0.10317987 1.91720424
31 32 33 34 35 36
-11.79254212 -1.44231029 -6.45741887 0.79536963 -4.51571225 -5.09551971
37 38 39 40 41 42
1.45965479 -3.30575920 -10.20577389 6.22426889 0.74971083 -9.21728670
43 44 45 46 47 48
0.88119446 -1.97949839 -1.76155409 -0.43351290 -4.37471135 -3.63442812
49 50 51 52 53 54
2.06734545 -7.13537855 5.02419156 1.97861342 -1.38076069 -5.48087221
55 56 57 58 59 60
1.19558753 4.45596130 2.70647158 3.50375984 3.43060694 6.24775772
61 62 63 64 65 66
2.62061391 -2.95562188 6.32384224 1.88910336 1.04766613 -6.04349616
67 68 69 70
0.11197908 3.48589698 5.83131926 -0.29754285
> postscript(file="/var/www/html/rcomp/tmp/6dlht1290976213.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 70
Frequency = 1
lag(myerror, k = 1) myerror
0 -4.13398365 NA
1 4.25096541 -4.13398365
2 2.17182653 4.25096541
3 -4.95742384 2.17182653
4 0.56718862 -4.95742384
5 10.89539276 0.56718862
6 11.62859205 10.89539276
7 -3.42459841 11.62859205
8 0.01524952 -3.42459841
9 -3.34739188 0.01524952
10 4.49170561 -3.34739188
11 2.85858576 4.49170561
12 -3.34154645 2.85858576
13 4.55507904 -3.34154645
14 -0.76785897 4.55507904
15 -2.77287537 -0.76785897
16 -0.88062501 -2.77287537
17 7.92905807 -0.88062501
18 -2.02481099 7.92905807
19 -1.09545118 -2.02481099
20 -0.33406741 -1.09545118
21 -0.22068185 -0.33406741
22 0.96811105 -0.22068185
23 -0.37639566 0.96811105
24 1.32791596 -0.37639566
25 4.59071519 1.32791596
26 -2.54622747 4.59071519
27 -2.36168647 -2.54622747
28 -0.10317987 -2.36168647
29 1.91720424 -0.10317987
30 -11.79254212 1.91720424
31 -1.44231029 -11.79254212
32 -6.45741887 -1.44231029
33 0.79536963 -6.45741887
34 -4.51571225 0.79536963
35 -5.09551971 -4.51571225
36 1.45965479 -5.09551971
37 -3.30575920 1.45965479
38 -10.20577389 -3.30575920
39 6.22426889 -10.20577389
40 0.74971083 6.22426889
41 -9.21728670 0.74971083
42 0.88119446 -9.21728670
43 -1.97949839 0.88119446
44 -1.76155409 -1.97949839
45 -0.43351290 -1.76155409
46 -4.37471135 -0.43351290
47 -3.63442812 -4.37471135
48 2.06734545 -3.63442812
49 -7.13537855 2.06734545
50 5.02419156 -7.13537855
51 1.97861342 5.02419156
52 -1.38076069 1.97861342
53 -5.48087221 -1.38076069
54 1.19558753 -5.48087221
55 4.45596130 1.19558753
56 2.70647158 4.45596130
57 3.50375984 2.70647158
58 3.43060694 3.50375984
59 6.24775772 3.43060694
60 2.62061391 6.24775772
61 -2.95562188 2.62061391
62 6.32384224 -2.95562188
63 1.88910336 6.32384224
64 1.04766613 1.88910336
65 -6.04349616 1.04766613
66 0.11197908 -6.04349616
67 3.48589698 0.11197908
68 5.83131926 3.48589698
69 -0.29754285 5.83131926
70 NA -0.29754285
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 4.25096541 -4.13398365
[2,] 2.17182653 4.25096541
[3,] -4.95742384 2.17182653
[4,] 0.56718862 -4.95742384
[5,] 10.89539276 0.56718862
[6,] 11.62859205 10.89539276
[7,] -3.42459841 11.62859205
[8,] 0.01524952 -3.42459841
[9,] -3.34739188 0.01524952
[10,] 4.49170561 -3.34739188
[11,] 2.85858576 4.49170561
[12,] -3.34154645 2.85858576
[13,] 4.55507904 -3.34154645
[14,] -0.76785897 4.55507904
[15,] -2.77287537 -0.76785897
[16,] -0.88062501 -2.77287537
[17,] 7.92905807 -0.88062501
[18,] -2.02481099 7.92905807
[19,] -1.09545118 -2.02481099
[20,] -0.33406741 -1.09545118
[21,] -0.22068185 -0.33406741
[22,] 0.96811105 -0.22068185
[23,] -0.37639566 0.96811105
[24,] 1.32791596 -0.37639566
[25,] 4.59071519 1.32791596
[26,] -2.54622747 4.59071519
[27,] -2.36168647 -2.54622747
[28,] -0.10317987 -2.36168647
[29,] 1.91720424 -0.10317987
[30,] -11.79254212 1.91720424
[31,] -1.44231029 -11.79254212
[32,] -6.45741887 -1.44231029
[33,] 0.79536963 -6.45741887
[34,] -4.51571225 0.79536963
[35,] -5.09551971 -4.51571225
[36,] 1.45965479 -5.09551971
[37,] -3.30575920 1.45965479
[38,] -10.20577389 -3.30575920
[39,] 6.22426889 -10.20577389
[40,] 0.74971083 6.22426889
[41,] -9.21728670 0.74971083
[42,] 0.88119446 -9.21728670
[43,] -1.97949839 0.88119446
[44,] -1.76155409 -1.97949839
[45,] -0.43351290 -1.76155409
[46,] -4.37471135 -0.43351290
[47,] -3.63442812 -4.37471135
[48,] 2.06734545 -3.63442812
[49,] -7.13537855 2.06734545
[50,] 5.02419156 -7.13537855
[51,] 1.97861342 5.02419156
[52,] -1.38076069 1.97861342
[53,] -5.48087221 -1.38076069
[54,] 1.19558753 -5.48087221
[55,] 4.45596130 1.19558753
[56,] 2.70647158 4.45596130
[57,] 3.50375984 2.70647158
[58,] 3.43060694 3.50375984
[59,] 6.24775772 3.43060694
[60,] 2.62061391 6.24775772
[61,] -2.95562188 2.62061391
[62,] 6.32384224 -2.95562188
[63,] 1.88910336 6.32384224
[64,] 1.04766613 1.88910336
[65,] -6.04349616 1.04766613
[66,] 0.11197908 -6.04349616
[67,] 3.48589698 0.11197908
[68,] 5.83131926 3.48589698
[69,] -0.29754285 5.83131926
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 4.25096541 -4.13398365
2 2.17182653 4.25096541
3 -4.95742384 2.17182653
4 0.56718862 -4.95742384
5 10.89539276 0.56718862
6 11.62859205 10.89539276
7 -3.42459841 11.62859205
8 0.01524952 -3.42459841
9 -3.34739188 0.01524952
10 4.49170561 -3.34739188
11 2.85858576 4.49170561
12 -3.34154645 2.85858576
13 4.55507904 -3.34154645
14 -0.76785897 4.55507904
15 -2.77287537 -0.76785897
16 -0.88062501 -2.77287537
17 7.92905807 -0.88062501
18 -2.02481099 7.92905807
19 -1.09545118 -2.02481099
20 -0.33406741 -1.09545118
21 -0.22068185 -0.33406741
22 0.96811105 -0.22068185
23 -0.37639566 0.96811105
24 1.32791596 -0.37639566
25 4.59071519 1.32791596
26 -2.54622747 4.59071519
27 -2.36168647 -2.54622747
28 -0.10317987 -2.36168647
29 1.91720424 -0.10317987
30 -11.79254212 1.91720424
31 -1.44231029 -11.79254212
32 -6.45741887 -1.44231029
33 0.79536963 -6.45741887
34 -4.51571225 0.79536963
35 -5.09551971 -4.51571225
36 1.45965479 -5.09551971
37 -3.30575920 1.45965479
38 -10.20577389 -3.30575920
39 6.22426889 -10.20577389
40 0.74971083 6.22426889
41 -9.21728670 0.74971083
42 0.88119446 -9.21728670
43 -1.97949839 0.88119446
44 -1.76155409 -1.97949839
45 -0.43351290 -1.76155409
46 -4.37471135 -0.43351290
47 -3.63442812 -4.37471135
48 2.06734545 -3.63442812
49 -7.13537855 2.06734545
50 5.02419156 -7.13537855
51 1.97861342 5.02419156
52 -1.38076069 1.97861342
53 -5.48087221 -1.38076069
54 1.19558753 -5.48087221
55 4.45596130 1.19558753
56 2.70647158 4.45596130
57 3.50375984 2.70647158
58 3.43060694 3.50375984
59 6.24775772 3.43060694
60 2.62061391 6.24775772
61 -2.95562188 2.62061391
62 6.32384224 -2.95562188
63 1.88910336 6.32384224
64 1.04766613 1.88910336
65 -6.04349616 1.04766613
66 0.11197908 -6.04349616
67 3.48589698 0.11197908
68 5.83131926 3.48589698
69 -0.29754285 5.83131926
> 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/76dyw1290976213.ps",horizontal=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/86dyw1290976213.ps",horizontal=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/96dyw1290976213.ps",horizontal=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/10gmfz1290976213.ps",horizontal=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/112me51290976213.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/1255us1290976213.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/13kfa11290976213.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/14nx8p1290976213.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/158g7v1290976213.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/16uyo11290976213.tab")
+ }
>
> try(system("convert tmp/1sli51290976213.ps tmp/1sli51290976213.png",intern=TRUE))
character(0)
> try(system("convert tmp/2sli51290976213.ps tmp/2sli51290976213.png",intern=TRUE))
character(0)
> try(system("convert tmp/32c0q1290976213.ps tmp/32c0q1290976213.png",intern=TRUE))
character(0)
> try(system("convert tmp/42c0q1290976213.ps tmp/42c0q1290976213.png",intern=TRUE))
character(0)
> try(system("convert tmp/52c0q1290976213.ps tmp/52c0q1290976213.png",intern=TRUE))
character(0)
> try(system("convert tmp/6dlht1290976213.ps tmp/6dlht1290976213.png",intern=TRUE))
character(0)
> try(system("convert tmp/76dyw1290976213.ps tmp/76dyw1290976213.png",intern=TRUE))
character(0)
> try(system("convert tmp/86dyw1290976213.ps tmp/86dyw1290976213.png",intern=TRUE))
character(0)
> try(system("convert tmp/96dyw1290976213.ps tmp/96dyw1290976213.png",intern=TRUE))
character(0)
> try(system("convert tmp/10gmfz1290976213.ps tmp/10gmfz1290976213.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.485 1.552 5.730