R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i686-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(102007
+ ,24776
+ ,4
+ ,112007
+ ,19814
+ ,4
+ ,122007
+ ,12738
+ ,4
+ ,012008
+ ,31566
+ ,4
+ ,022008
+ ,30111
+ ,4
+ ,032008
+ ,30019
+ ,4
+ ,042008
+ ,31934
+ ,4
+ ,052008
+ ,25826
+ ,4
+ ,062008
+ ,26835
+ ,4
+ ,072008
+ ,20205
+ ,4.18
+ ,082008
+ ,17789
+ ,4.25
+ ,092008
+ ,20520
+ ,4.25
+ ,102008
+ ,22518
+ ,3.97
+ ,112008
+ ,15572
+ ,3.42
+ ,122008
+ ,11509
+ ,2.75
+ ,012009
+ ,25447
+ ,2.31
+ ,022009
+ ,24090
+ ,2
+ ,032009
+ ,27786
+ ,1.66
+ ,042009
+ ,26195
+ ,1.31
+ ,052009
+ ,20516
+ ,1.09
+ ,062009
+ ,22759
+ ,1
+ ,072009
+ ,19028
+ ,1
+ ,082009
+ ,16971
+ ,1
+ ,092009
+ ,20036
+ ,1
+ ,102009
+ ,22485
+ ,1
+ ,112009
+ ,18730
+ ,1
+ ,122009
+ ,14538
+ ,1
+ ,012010
+ ,27561
+ ,1
+ ,022010
+ ,25985
+ ,1
+ ,032010
+ ,34670
+ ,1
+ ,042010
+ ,32066
+ ,1
+ ,052010
+ ,27186
+ ,1
+ ,062010
+ ,29586
+ ,1
+ ,072010
+ ,21359
+ ,1
+ ,082010
+ ,21553
+ ,1
+ ,092010
+ ,19573
+ ,1
+ ,102010
+ ,24256
+ ,1
+ ,112010
+ ,22380
+ ,1
+ ,122010
+ ,16167
+ ,1
+ ,012011
+ ,27297
+ ,1
+ ,022011
+ ,28287
+ ,1
+ ,032011
+ ,33474
+ ,1
+ ,042011
+ ,28229
+ ,1.14
+ ,052011
+ ,28785
+ ,1.25
+ ,062011
+ ,25597
+ ,1.25
+ ,072011
+ ,18130
+ ,1.4
+ ,082011
+ ,20198
+ ,1.5
+ ,092011
+ ,22849
+ ,1.5
+ ,102011
+ ,23118
+ ,1.5
+ ,112011
+ ,21925
+ ,1.32
+ ,122011
+ ,20801
+ ,1.11
+ ,012012
+ ,18785
+ ,1
+ ,122012
+ ,20659
+ ,1
+ ,032012
+ ,29367
+ ,1
+ ,042012
+ ,23992
+ ,1
+ ,052012
+ ,20645
+ ,1
+ ,062012
+ ,22356
+ ,1
+ ,072012
+ ,17902
+ ,0.83
+ ,082012
+ ,15879
+ ,0.75
+ ,092012
+ ,16963
+ ,0.75
+ ,102012
+ ,21035
+ ,0.75)
+ ,dim=c(3
+ ,61)
+ ,dimnames=list(c('Data'
+ ,'Inschrijvingen'
+ ,'Rentevoet')
+ ,1:61))
> y <- array(NA,dim=c(3,61),dimnames=list(c('Data','Inschrijvingen','Rentevoet'),1:61))
> 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 = 'Do not include Seasonal Dummies'
> par1 = '2'
> par3 <- 'No Linear Trend'
> par2 <- 'Do not include Seasonal Dummies'
> par1 <- '2'
> #'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, 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
Inschrijvingen Data Rentevoet
1 24776 102007 4.00
2 19814 112007 4.00
3 12738 122007 4.00
4 31566 12008 4.00
5 30111 22008 4.00
6 30019 32008 4.00
7 31934 42008 4.00
8 25826 52008 4.00
9 26835 62008 4.00
10 20205 72008 4.18
11 17789 82008 4.25
12 20520 92008 4.25
13 22518 102008 3.97
14 15572 112008 3.42
15 11509 122008 2.75
16 25447 12009 2.31
17 24090 22009 2.00
18 27786 32009 1.66
19 26195 42009 1.31
20 20516 52009 1.09
21 22759 62009 1.00
22 19028 72009 1.00
23 16971 82009 1.00
24 20036 92009 1.00
25 22485 102009 1.00
26 18730 112009 1.00
27 14538 122009 1.00
28 27561 12010 1.00
29 25985 22010 1.00
30 34670 32010 1.00
31 32066 42010 1.00
32 27186 52010 1.00
33 29586 62010 1.00
34 21359 72010 1.00
35 21553 82010 1.00
36 19573 92010 1.00
37 24256 102010 1.00
38 22380 112010 1.00
39 16167 122010 1.00
40 27297 12011 1.00
41 28287 22011 1.00
42 33474 32011 1.00
43 28229 42011 1.14
44 28785 52011 1.25
45 25597 62011 1.25
46 18130 72011 1.40
47 20198 82011 1.50
48 22849 92011 1.50
49 23118 102011 1.50
50 21925 112011 1.32
51 20801 122011 1.11
52 18785 12012 1.00
53 20659 122012 1.00
54 29367 32012 1.00
55 23992 42012 1.00
56 20645 52012 1.00
57 22356 62012 1.00
58 17902 72012 0.83
59 15879 82012 0.75
60 16963 92012 0.75
61 21035 102012 0.75
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Data Rentevoet
29959.2732 -0.1055 245.4704
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-10152.0 -3180.7 346.6 2470.3 7843.6
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.996e+04 1.266e+03 23.657 < 2e-16 ***
Data -1.055e-01 1.392e-02 -7.581 3.09e-10 ***
Rentevoet 2.455e+02 3.873e+02 0.634 0.529
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3770 on 58 degrees of freedom
Multiple R-squared: 0.4983, Adjusted R-squared: 0.481
F-statistic: 28.8 on 2 and 58 DF, p-value: 2.056e-09
> 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.62390475 0.75219049 0.37609525
[2,] 0.60061869 0.79876261 0.39938131
[3,] 0.47778916 0.95557832 0.52221084
[4,] 0.36909033 0.73818065 0.63090967
[5,] 0.25496344 0.50992687 0.74503656
[6,] 0.17049027 0.34098053 0.82950973
[7,] 0.16666299 0.33332599 0.83333701
[8,] 0.13056621 0.26113242 0.86943379
[9,] 0.27305566 0.54611132 0.72694434
[10,] 0.26986091 0.53972182 0.73013909
[11,] 0.21870612 0.43741224 0.78129388
[12,] 0.18229411 0.36458822 0.81770589
[13,] 0.23043375 0.46086750 0.76956625
[14,] 0.22307968 0.44615936 0.77692032
[15,] 0.18526248 0.37052496 0.81473752
[16,] 0.15871541 0.31743082 0.84128459
[17,] 0.12473116 0.24946231 0.87526884
[18,] 0.10616334 0.21232668 0.89383666
[19,] 0.10155852 0.20311704 0.89844148
[20,] 0.16919485 0.33838970 0.83080515
[21,] 0.14619003 0.29238007 0.85380997
[22,] 0.11957123 0.23914245 0.88042877
[23,] 0.08577151 0.17154303 0.91422849
[24,] 0.06155886 0.12311773 0.93844114
[25,] 0.25046066 0.50092133 0.74953934
[26,] 0.41146143 0.82292287 0.58853857
[27,] 0.38192276 0.76384553 0.61807724
[28,] 0.52552252 0.94895496 0.47447748
[29,] 0.45303420 0.90606841 0.54696580
[30,] 0.37787224 0.75574448 0.62212776
[31,] 0.31106170 0.62212340 0.68893830
[32,] 0.36168757 0.72337513 0.63831243
[33,] 0.36934822 0.73869645 0.63065178
[34,] 0.30759329 0.61518658 0.69240671
[35,] 0.25336525 0.50673051 0.74663475
[36,] 0.20755604 0.41511207 0.79244396
[37,] 0.50293429 0.99413142 0.49706571
[38,] 0.53093349 0.93813303 0.46906651
[39,] 0.64219661 0.71560679 0.35780339
[40,] 0.62713628 0.74572744 0.37286372
[41,] 0.70138644 0.59722711 0.29861356
[42,] 0.70205970 0.59588061 0.29794030
[43,] 0.63093425 0.73813150 0.36906575
[44,] 0.56792362 0.86415277 0.43207638
[45,] 0.51188989 0.97622022 0.48811011
[46,] 0.44153951 0.88307903 0.55846049
[47,] 0.68722879 0.62554242 0.31277121
[48,] 0.57089128 0.85821745 0.42910872
[49,] 0.86748140 0.26503720 0.13251860
[50,] 0.98237602 0.03524796 0.01762398
> postscript(file="/var/wessaorg/rcomp/tmp/1k09g1353332705.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/2ktp11353332705.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/32db81353332705.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/44qsv1353332705.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/56wvl1353332705.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 = 61
Frequency = 1
1 2 3 4 5
4600.553728 693.942901 -5326.667926 1892.156707 1492.545880
6 7 8 9 10
2455.935054 5426.324227 373.713400 2438.102573 -3180.692919
11 12 13 14 15
-4558.486671 -772.097497 2350.023377 -3405.578750 -6248.724434
16 17 18 19 20
-3811.892842 -4037.407856 797.441240 347.745040 -4221.862307
21 22 23 24 25
-901.380801 -3576.991628 -4578.602454 -458.213281 3046.175892
26 27 28 29 30
346.565065 -2790.045761 -1376.221128 -1896.831955 7843.557218
31 32 33 34 35
6294.946392 2470.335565 5925.724738 -1245.886089 3.503085
36 37 38 39 40
-921.107742 4817.281431 3996.670604 -1160.940222 -1640.115589
41 42 43 44 45
405.273584 6647.662757 2423.686080 4008.073513 1875.462686
46 47 48 49 50
-4572.968695 -1474.126558 2232.262616 3556.651789 3463.225627
51 52 53 54 55
3446.163577 -10152.010050 3331.270855 2540.768296 -1778.842531
56 57 58 59 60
-4070.453357 -1304.064184 -4660.945049 -5608.918247 -3469.529074
61
1657.860099
> postscript(file="/var/wessaorg/rcomp/tmp/6d6o11353332705.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 = 61
Frequency = 1
lag(myerror, k = 1) myerror
0 4600.553728 NA
1 693.942901 4600.553728
2 -5326.667926 693.942901
3 1892.156707 -5326.667926
4 1492.545880 1892.156707
5 2455.935054 1492.545880
6 5426.324227 2455.935054
7 373.713400 5426.324227
8 2438.102573 373.713400
9 -3180.692919 2438.102573
10 -4558.486671 -3180.692919
11 -772.097497 -4558.486671
12 2350.023377 -772.097497
13 -3405.578750 2350.023377
14 -6248.724434 -3405.578750
15 -3811.892842 -6248.724434
16 -4037.407856 -3811.892842
17 797.441240 -4037.407856
18 347.745040 797.441240
19 -4221.862307 347.745040
20 -901.380801 -4221.862307
21 -3576.991628 -901.380801
22 -4578.602454 -3576.991628
23 -458.213281 -4578.602454
24 3046.175892 -458.213281
25 346.565065 3046.175892
26 -2790.045761 346.565065
27 -1376.221128 -2790.045761
28 -1896.831955 -1376.221128
29 7843.557218 -1896.831955
30 6294.946392 7843.557218
31 2470.335565 6294.946392
32 5925.724738 2470.335565
33 -1245.886089 5925.724738
34 3.503085 -1245.886089
35 -921.107742 3.503085
36 4817.281431 -921.107742
37 3996.670604 4817.281431
38 -1160.940222 3996.670604
39 -1640.115589 -1160.940222
40 405.273584 -1640.115589
41 6647.662757 405.273584
42 2423.686080 6647.662757
43 4008.073513 2423.686080
44 1875.462686 4008.073513
45 -4572.968695 1875.462686
46 -1474.126558 -4572.968695
47 2232.262616 -1474.126558
48 3556.651789 2232.262616
49 3463.225627 3556.651789
50 3446.163577 3463.225627
51 -10152.010050 3446.163577
52 3331.270855 -10152.010050
53 2540.768296 3331.270855
54 -1778.842531 2540.768296
55 -4070.453357 -1778.842531
56 -1304.064184 -4070.453357
57 -4660.945049 -1304.064184
58 -5608.918247 -4660.945049
59 -3469.529074 -5608.918247
60 1657.860099 -3469.529074
61 NA 1657.860099
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 693.942901 4600.553728
[2,] -5326.667926 693.942901
[3,] 1892.156707 -5326.667926
[4,] 1492.545880 1892.156707
[5,] 2455.935054 1492.545880
[6,] 5426.324227 2455.935054
[7,] 373.713400 5426.324227
[8,] 2438.102573 373.713400
[9,] -3180.692919 2438.102573
[10,] -4558.486671 -3180.692919
[11,] -772.097497 -4558.486671
[12,] 2350.023377 -772.097497
[13,] -3405.578750 2350.023377
[14,] -6248.724434 -3405.578750
[15,] -3811.892842 -6248.724434
[16,] -4037.407856 -3811.892842
[17,] 797.441240 -4037.407856
[18,] 347.745040 797.441240
[19,] -4221.862307 347.745040
[20,] -901.380801 -4221.862307
[21,] -3576.991628 -901.380801
[22,] -4578.602454 -3576.991628
[23,] -458.213281 -4578.602454
[24,] 3046.175892 -458.213281
[25,] 346.565065 3046.175892
[26,] -2790.045761 346.565065
[27,] -1376.221128 -2790.045761
[28,] -1896.831955 -1376.221128
[29,] 7843.557218 -1896.831955
[30,] 6294.946392 7843.557218
[31,] 2470.335565 6294.946392
[32,] 5925.724738 2470.335565
[33,] -1245.886089 5925.724738
[34,] 3.503085 -1245.886089
[35,] -921.107742 3.503085
[36,] 4817.281431 -921.107742
[37,] 3996.670604 4817.281431
[38,] -1160.940222 3996.670604
[39,] -1640.115589 -1160.940222
[40,] 405.273584 -1640.115589
[41,] 6647.662757 405.273584
[42,] 2423.686080 6647.662757
[43,] 4008.073513 2423.686080
[44,] 1875.462686 4008.073513
[45,] -4572.968695 1875.462686
[46,] -1474.126558 -4572.968695
[47,] 2232.262616 -1474.126558
[48,] 3556.651789 2232.262616
[49,] 3463.225627 3556.651789
[50,] 3446.163577 3463.225627
[51,] -10152.010050 3446.163577
[52,] 3331.270855 -10152.010050
[53,] 2540.768296 3331.270855
[54,] -1778.842531 2540.768296
[55,] -4070.453357 -1778.842531
[56,] -1304.064184 -4070.453357
[57,] -4660.945049 -1304.064184
[58,] -5608.918247 -4660.945049
[59,] -3469.529074 -5608.918247
[60,] 1657.860099 -3469.529074
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 693.942901 4600.553728
2 -5326.667926 693.942901
3 1892.156707 -5326.667926
4 1492.545880 1892.156707
5 2455.935054 1492.545880
6 5426.324227 2455.935054
7 373.713400 5426.324227
8 2438.102573 373.713400
9 -3180.692919 2438.102573
10 -4558.486671 -3180.692919
11 -772.097497 -4558.486671
12 2350.023377 -772.097497
13 -3405.578750 2350.023377
14 -6248.724434 -3405.578750
15 -3811.892842 -6248.724434
16 -4037.407856 -3811.892842
17 797.441240 -4037.407856
18 347.745040 797.441240
19 -4221.862307 347.745040
20 -901.380801 -4221.862307
21 -3576.991628 -901.380801
22 -4578.602454 -3576.991628
23 -458.213281 -4578.602454
24 3046.175892 -458.213281
25 346.565065 3046.175892
26 -2790.045761 346.565065
27 -1376.221128 -2790.045761
28 -1896.831955 -1376.221128
29 7843.557218 -1896.831955
30 6294.946392 7843.557218
31 2470.335565 6294.946392
32 5925.724738 2470.335565
33 -1245.886089 5925.724738
34 3.503085 -1245.886089
35 -921.107742 3.503085
36 4817.281431 -921.107742
37 3996.670604 4817.281431
38 -1160.940222 3996.670604
39 -1640.115589 -1160.940222
40 405.273584 -1640.115589
41 6647.662757 405.273584
42 2423.686080 6647.662757
43 4008.073513 2423.686080
44 1875.462686 4008.073513
45 -4572.968695 1875.462686
46 -1474.126558 -4572.968695
47 2232.262616 -1474.126558
48 3556.651789 2232.262616
49 3463.225627 3556.651789
50 3446.163577 3463.225627
51 -10152.010050 3446.163577
52 3331.270855 -10152.010050
53 2540.768296 3331.270855
54 -1778.842531 2540.768296
55 -4070.453357 -1778.842531
56 -1304.064184 -4070.453357
57 -4660.945049 -1304.064184
58 -5608.918247 -4660.945049
59 -3469.529074 -5608.918247
60 1657.860099 -3469.529074
> 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/76cp51353332705.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/83p761353332705.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/9ldri1353332705.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/10jvcv1353332705.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/11evpj1353332705.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/12xbwc1353332705.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/13nu6l1353332705.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/14c8f31353332705.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/1535kv1353332705.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/16iin91353332706.tab")
+ }
>
> try(system("convert tmp/1k09g1353332705.ps tmp/1k09g1353332705.png",intern=TRUE))
character(0)
> try(system("convert tmp/2ktp11353332705.ps tmp/2ktp11353332705.png",intern=TRUE))
character(0)
> try(system("convert tmp/32db81353332705.ps tmp/32db81353332705.png",intern=TRUE))
character(0)
> try(system("convert tmp/44qsv1353332705.ps tmp/44qsv1353332705.png",intern=TRUE))
character(0)
> try(system("convert tmp/56wvl1353332705.ps tmp/56wvl1353332705.png",intern=TRUE))
character(0)
> try(system("convert tmp/6d6o11353332705.ps tmp/6d6o11353332705.png",intern=TRUE))
character(0)
> try(system("convert tmp/76cp51353332705.ps tmp/76cp51353332705.png",intern=TRUE))
character(0)
> try(system("convert tmp/83p761353332705.ps tmp/83p761353332705.png",intern=TRUE))
character(0)
> try(system("convert tmp/9ldri1353332705.ps tmp/9ldri1353332705.png",intern=TRUE))
character(0)
> try(system("convert tmp/10jvcv1353332705.ps tmp/10jvcv1353332705.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
9.735 1.806 15.156