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(4.1,96.8,96.2,3.9,109.9,96.8,3.8,88,109.9,3.7,91.1,88,3.7,106.4,91.1,4.1,68.6,106.4,4.1,100.1,68.6,3.8,108,100.1,3.7,106,108,3.5,108.6,106,3.6,91.5,108.6,4.1,99.2,91.5,3.8,98,99.2,3.7,96.6,98,3.6,102.8,96.6,3.3,96.9,102.8,3.4,110,96.9,3.7,70.5,110,3.7,101.9,70.5,3.4,109.6,101.9,3.3,107.8,109.6,3,113,107.8,3,93.8,113,3.3,108,93.8,3,102.8,108,2.9,116.3,102.8,2.8,89.2,116.3,2.5,106.7,89.2,2.6,112.1,106.7,2.8,74.2,112.1,2.7,108.8,74.2,2.4,111.5,108.8,2.2,118.8,111.5,2.1,118.9,118.8,2.1,97.6,118.9,2.3,116.4,97.6,2.1,107.9,116.4,2,121.2,107.9,1.9,97.9,121.2,1.7,113.4,97.9,1.8,117.6,113.4,2.1,79.6,117.6,2,115.9,79.6,1.8,115.7,115.9,1.7,129.1,115.7,1.6,123.3,129.1,1.6,96.7,123.3,1.8,121.2,96.7,1.7,118.2,121.2,1.7,102.1,118.2,1.5,125.4,102.1,1.5,116.7,125.4,1.5,121.3,116.7,1.8,85.3,121.3,1.8,114.2,85.3,1.7,124.4,114.2,1.7,131,124.4,1.8,118.3,131,2,99.6,118.3),dim=c(3,59),dimnames=list(c('unempl','proman','Y(t-1)'),1:59))
> y <- array(NA,dim=c(3,59),dimnames=list(c('unempl','proman','Y(t-1)'),1:59))
> 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
unempl proman Y(t-1) M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 4.1 96.8 96.2 1 0 0 0 0 0 0 0 0 0 0 1
2 3.9 109.9 96.8 0 1 0 0 0 0 0 0 0 0 0 2
3 3.8 88.0 109.9 0 0 1 0 0 0 0 0 0 0 0 3
4 3.7 91.1 88.0 0 0 0 1 0 0 0 0 0 0 0 4
5 3.7 106.4 91.1 0 0 0 0 1 0 0 0 0 0 0 5
6 4.1 68.6 106.4 0 0 0 0 0 1 0 0 0 0 0 6
7 4.1 100.1 68.6 0 0 0 0 0 0 1 0 0 0 0 7
8 3.8 108.0 100.1 0 0 0 0 0 0 0 1 0 0 0 8
9 3.7 106.0 108.0 0 0 0 0 0 0 0 0 1 0 0 9
10 3.5 108.6 106.0 0 0 0 0 0 0 0 0 0 1 0 10
11 3.6 91.5 108.6 0 0 0 0 0 0 0 0 0 0 1 11
12 4.1 99.2 91.5 0 0 0 0 0 0 0 0 0 0 0 12
13 3.8 98.0 99.2 1 0 0 0 0 0 0 0 0 0 0 13
14 3.7 96.6 98.0 0 1 0 0 0 0 0 0 0 0 0 14
15 3.6 102.8 96.6 0 0 1 0 0 0 0 0 0 0 0 15
16 3.3 96.9 102.8 0 0 0 1 0 0 0 0 0 0 0 16
17 3.4 110.0 96.9 0 0 0 0 1 0 0 0 0 0 0 17
18 3.7 70.5 110.0 0 0 0 0 0 1 0 0 0 0 0 18
19 3.7 101.9 70.5 0 0 0 0 0 0 1 0 0 0 0 19
20 3.4 109.6 101.9 0 0 0 0 0 0 0 1 0 0 0 20
21 3.3 107.8 109.6 0 0 0 0 0 0 0 0 1 0 0 21
22 3.0 113.0 107.8 0 0 0 0 0 0 0 0 0 1 0 22
23 3.0 93.8 113.0 0 0 0 0 0 0 0 0 0 0 1 23
24 3.3 108.0 93.8 0 0 0 0 0 0 0 0 0 0 0 24
25 3.0 102.8 108.0 1 0 0 0 0 0 0 0 0 0 0 25
26 2.9 116.3 102.8 0 1 0 0 0 0 0 0 0 0 0 26
27 2.8 89.2 116.3 0 0 1 0 0 0 0 0 0 0 0 27
28 2.5 106.7 89.2 0 0 0 1 0 0 0 0 0 0 0 28
29 2.6 112.1 106.7 0 0 0 0 1 0 0 0 0 0 0 29
30 2.8 74.2 112.1 0 0 0 0 0 1 0 0 0 0 0 30
31 2.7 108.8 74.2 0 0 0 0 0 0 1 0 0 0 0 31
32 2.4 111.5 108.8 0 0 0 0 0 0 0 1 0 0 0 32
33 2.2 118.8 111.5 0 0 0 0 0 0 0 0 1 0 0 33
34 2.1 118.9 118.8 0 0 0 0 0 0 0 0 0 1 0 34
35 2.1 97.6 118.9 0 0 0 0 0 0 0 0 0 0 1 35
36 2.3 116.4 97.6 0 0 0 0 0 0 0 0 0 0 0 36
37 2.1 107.9 116.4 1 0 0 0 0 0 0 0 0 0 0 37
38 2.0 121.2 107.9 0 1 0 0 0 0 0 0 0 0 0 38
39 1.9 97.9 121.2 0 0 1 0 0 0 0 0 0 0 0 39
40 1.7 113.4 97.9 0 0 0 1 0 0 0 0 0 0 0 40
41 1.8 117.6 113.4 0 0 0 0 1 0 0 0 0 0 0 41
42 2.1 79.6 117.6 0 0 0 0 0 1 0 0 0 0 0 42
43 2.0 115.9 79.6 0 0 0 0 0 0 1 0 0 0 0 43
44 1.8 115.7 115.9 0 0 0 0 0 0 0 1 0 0 0 44
45 1.7 129.1 115.7 0 0 0 0 0 0 0 0 1 0 0 45
46 1.6 123.3 129.1 0 0 0 0 0 0 0 0 0 1 0 46
47 1.6 96.7 123.3 0 0 0 0 0 0 0 0 0 0 1 47
48 1.8 121.2 96.7 0 0 0 0 0 0 0 0 0 0 0 48
49 1.7 118.2 121.2 1 0 0 0 0 0 0 0 0 0 0 49
50 1.7 102.1 118.2 0 1 0 0 0 0 0 0 0 0 0 50
51 1.5 125.4 102.1 0 0 1 0 0 0 0 0 0 0 0 51
52 1.5 116.7 125.4 0 0 0 1 0 0 0 0 0 0 0 52
53 1.5 121.3 116.7 0 0 0 0 1 0 0 0 0 0 0 53
54 1.8 85.3 121.3 0 0 0 0 0 1 0 0 0 0 0 54
55 1.8 114.2 85.3 0 0 0 0 0 0 1 0 0 0 0 55
56 1.7 124.4 114.2 0 0 0 0 0 0 0 1 0 0 0 56
57 1.7 131.0 124.4 0 0 0 0 0 0 0 0 1 0 0 57
58 1.8 118.3 131.0 0 0 0 0 0 0 0 0 0 1 0 58
59 2.0 99.6 118.3 0 0 0 0 0 0 0 0 0 0 1 59
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) proman `Y(t-1)` M1 M2 M3
6.958634 -0.017226 -0.010582 -0.099512 -0.120160 -0.301406
M4 M5 M6 M7 M8 M9
-0.459125 -0.168063 -0.390920 -0.232042 0.008602 0.088254
M10 M11 t
0.020269 -0.257871 -0.038796
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-0.340284 -0.169361 -0.005358 0.115598 0.555769
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.958634 1.085142 6.413 8.34e-08 ***
proman -0.017226 0.006944 -2.481 0.017 *
`Y(t-1)` -0.010582 0.007009 -1.510 0.138
M1 -0.099512 0.189178 -0.526 0.602
M2 -0.120160 0.179817 -0.668 0.507
M3 -0.301406 0.190444 -1.583 0.121
M4 -0.459125 0.167970 -2.733 0.009 **
M5 -0.168063 0.181497 -0.926 0.360
M6 -0.390920 0.277812 -1.407 0.166
M7 -0.232042 0.219876 -1.055 0.297
M8 0.008602 0.186889 0.046 0.963
M9 0.088254 0.219093 0.403 0.689
M10 0.020269 0.232632 0.087 0.931
M11 -0.257871 0.213585 -1.207 0.234
t -0.038796 0.004869 -7.968 4.48e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.2409 on 44 degrees of freedom
Multiple R-squared: 0.9436, Adjusted R-squared: 0.9256
F-statistic: 52.57 on 14 and 44 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.0055146266 0.0110292531 0.9944854
[2,] 0.0029067489 0.0058134978 0.9970933
[3,] 0.0012297908 0.0024595816 0.9987702
[4,] 0.0004351577 0.0008703154 0.9995648
[5,] 0.0004574807 0.0009149614 0.9995425
[6,] 0.0010962012 0.0021924024 0.9989038
[7,] 0.0119834232 0.0239668465 0.9880166
[8,] 0.0153828749 0.0307657499 0.9846171
[9,] 0.0224352687 0.0448705374 0.9775647
[10,] 0.0301826301 0.0603652602 0.9698174
[11,] 0.0566120746 0.1132241492 0.9433879
[12,] 0.0673630801 0.1347261602 0.9326369
[13,] 0.1496172346 0.2992344692 0.8503828
[14,] 0.3040914196 0.6081828392 0.6959086
[15,] 0.4083361651 0.8166723302 0.5916638
[16,] 0.4260236790 0.8520473579 0.5739763
[17,] 0.3415667549 0.6831335099 0.6584332
[18,] 0.2917489755 0.5834979511 0.7082510
[19,] 0.3824814197 0.7649628394 0.6175186
[20,] 0.3166723818 0.6333447636 0.6833276
[21,] 0.5346340061 0.9307319878 0.4653660
[22,] 0.4085225438 0.8170450877 0.5914775
[23,] 0.5377819962 0.9244360077 0.4622180
[24,] 0.4695717100 0.9391434200 0.5304283
> postscript(file="/var/www/html/rcomp/tmp/1t1kx1258666076.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/2a7zg1258666076.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/39onr1258666076.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/49vfk1258666076.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/5fmmj1258666076.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 = 59
Frequency = 1
1 2 3 4 5 6
-0.034865824 0.056583030 -0.061985885 -0.143825548 -0.099735068 0.072702508
7 8 9 10 11 12
0.095211215 0.062790805 -0.028915951 -0.098513242 0.051380574 0.283983427
13 14 15 16 17 18
0.183105122 0.105734978 0.317759925 0.178255117 0.189207908 0.209080780
19 20 21 22 23 24
0.211876879 0.174953118 0.084574999 -0.038119344 0.003114978 0.125461008
25 26 27 28 29 30
0.024465925 0.161427186 -0.042481719 -0.131301619 -0.005357936 -0.139408503
31 32 33 34 35 36
-0.164558780 -0.253746748 -0.340284128 -0.254529177 -0.303438695 -0.224078064
37 38 39 40 41 42
-0.233238555 -0.134644270 -0.275212501 -0.258270479 -0.174162240 -0.222634223
43 44 45 46 47 48
-0.219559276 -0.240711335 -0.152861671 -0.104185026 -0.306826102 -0.185366372
49 50 51 52 53 54
0.060533331 -0.189100924 0.061920180 0.355142529 0.090047336 0.080259438
55 56 57 58 59
0.077029961 0.256714159 0.437486752 0.495346789 0.555769246
> postscript(file="/var/www/html/rcomp/tmp/6kvie1258666076.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 = 59
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.034865824 NA
1 0.056583030 -0.034865824
2 -0.061985885 0.056583030
3 -0.143825548 -0.061985885
4 -0.099735068 -0.143825548
5 0.072702508 -0.099735068
6 0.095211215 0.072702508
7 0.062790805 0.095211215
8 -0.028915951 0.062790805
9 -0.098513242 -0.028915951
10 0.051380574 -0.098513242
11 0.283983427 0.051380574
12 0.183105122 0.283983427
13 0.105734978 0.183105122
14 0.317759925 0.105734978
15 0.178255117 0.317759925
16 0.189207908 0.178255117
17 0.209080780 0.189207908
18 0.211876879 0.209080780
19 0.174953118 0.211876879
20 0.084574999 0.174953118
21 -0.038119344 0.084574999
22 0.003114978 -0.038119344
23 0.125461008 0.003114978
24 0.024465925 0.125461008
25 0.161427186 0.024465925
26 -0.042481719 0.161427186
27 -0.131301619 -0.042481719
28 -0.005357936 -0.131301619
29 -0.139408503 -0.005357936
30 -0.164558780 -0.139408503
31 -0.253746748 -0.164558780
32 -0.340284128 -0.253746748
33 -0.254529177 -0.340284128
34 -0.303438695 -0.254529177
35 -0.224078064 -0.303438695
36 -0.233238555 -0.224078064
37 -0.134644270 -0.233238555
38 -0.275212501 -0.134644270
39 -0.258270479 -0.275212501
40 -0.174162240 -0.258270479
41 -0.222634223 -0.174162240
42 -0.219559276 -0.222634223
43 -0.240711335 -0.219559276
44 -0.152861671 -0.240711335
45 -0.104185026 -0.152861671
46 -0.306826102 -0.104185026
47 -0.185366372 -0.306826102
48 0.060533331 -0.185366372
49 -0.189100924 0.060533331
50 0.061920180 -0.189100924
51 0.355142529 0.061920180
52 0.090047336 0.355142529
53 0.080259438 0.090047336
54 0.077029961 0.080259438
55 0.256714159 0.077029961
56 0.437486752 0.256714159
57 0.495346789 0.437486752
58 0.555769246 0.495346789
59 NA 0.555769246
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.056583030 -0.034865824
[2,] -0.061985885 0.056583030
[3,] -0.143825548 -0.061985885
[4,] -0.099735068 -0.143825548
[5,] 0.072702508 -0.099735068
[6,] 0.095211215 0.072702508
[7,] 0.062790805 0.095211215
[8,] -0.028915951 0.062790805
[9,] -0.098513242 -0.028915951
[10,] 0.051380574 -0.098513242
[11,] 0.283983427 0.051380574
[12,] 0.183105122 0.283983427
[13,] 0.105734978 0.183105122
[14,] 0.317759925 0.105734978
[15,] 0.178255117 0.317759925
[16,] 0.189207908 0.178255117
[17,] 0.209080780 0.189207908
[18,] 0.211876879 0.209080780
[19,] 0.174953118 0.211876879
[20,] 0.084574999 0.174953118
[21,] -0.038119344 0.084574999
[22,] 0.003114978 -0.038119344
[23,] 0.125461008 0.003114978
[24,] 0.024465925 0.125461008
[25,] 0.161427186 0.024465925
[26,] -0.042481719 0.161427186
[27,] -0.131301619 -0.042481719
[28,] -0.005357936 -0.131301619
[29,] -0.139408503 -0.005357936
[30,] -0.164558780 -0.139408503
[31,] -0.253746748 -0.164558780
[32,] -0.340284128 -0.253746748
[33,] -0.254529177 -0.340284128
[34,] -0.303438695 -0.254529177
[35,] -0.224078064 -0.303438695
[36,] -0.233238555 -0.224078064
[37,] -0.134644270 -0.233238555
[38,] -0.275212501 -0.134644270
[39,] -0.258270479 -0.275212501
[40,] -0.174162240 -0.258270479
[41,] -0.222634223 -0.174162240
[42,] -0.219559276 -0.222634223
[43,] -0.240711335 -0.219559276
[44,] -0.152861671 -0.240711335
[45,] -0.104185026 -0.152861671
[46,] -0.306826102 -0.104185026
[47,] -0.185366372 -0.306826102
[48,] 0.060533331 -0.185366372
[49,] -0.189100924 0.060533331
[50,] 0.061920180 -0.189100924
[51,] 0.355142529 0.061920180
[52,] 0.090047336 0.355142529
[53,] 0.080259438 0.090047336
[54,] 0.077029961 0.080259438
[55,] 0.256714159 0.077029961
[56,] 0.437486752 0.256714159
[57,] 0.495346789 0.437486752
[58,] 0.555769246 0.495346789
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.056583030 -0.034865824
2 -0.061985885 0.056583030
3 -0.143825548 -0.061985885
4 -0.099735068 -0.143825548
5 0.072702508 -0.099735068
6 0.095211215 0.072702508
7 0.062790805 0.095211215
8 -0.028915951 0.062790805
9 -0.098513242 -0.028915951
10 0.051380574 -0.098513242
11 0.283983427 0.051380574
12 0.183105122 0.283983427
13 0.105734978 0.183105122
14 0.317759925 0.105734978
15 0.178255117 0.317759925
16 0.189207908 0.178255117
17 0.209080780 0.189207908
18 0.211876879 0.209080780
19 0.174953118 0.211876879
20 0.084574999 0.174953118
21 -0.038119344 0.084574999
22 0.003114978 -0.038119344
23 0.125461008 0.003114978
24 0.024465925 0.125461008
25 0.161427186 0.024465925
26 -0.042481719 0.161427186
27 -0.131301619 -0.042481719
28 -0.005357936 -0.131301619
29 -0.139408503 -0.005357936
30 -0.164558780 -0.139408503
31 -0.253746748 -0.164558780
32 -0.340284128 -0.253746748
33 -0.254529177 -0.340284128
34 -0.303438695 -0.254529177
35 -0.224078064 -0.303438695
36 -0.233238555 -0.224078064
37 -0.134644270 -0.233238555
38 -0.275212501 -0.134644270
39 -0.258270479 -0.275212501
40 -0.174162240 -0.258270479
41 -0.222634223 -0.174162240
42 -0.219559276 -0.222634223
43 -0.240711335 -0.219559276
44 -0.152861671 -0.240711335
45 -0.104185026 -0.152861671
46 -0.306826102 -0.104185026
47 -0.185366372 -0.306826102
48 0.060533331 -0.185366372
49 -0.189100924 0.060533331
50 0.061920180 -0.189100924
51 0.355142529 0.061920180
52 0.090047336 0.355142529
53 0.080259438 0.090047336
54 0.077029961 0.080259438
55 0.256714159 0.077029961
56 0.437486752 0.256714159
57 0.495346789 0.437486752
58 0.555769246 0.495346789
> 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/7w07u1258666076.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/8wbxj1258666076.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/9myqd1258666076.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/10xqmf1258666076.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/11js5b1258666076.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/1272l31258666076.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/135ggn1258666076.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/14n6b41258666076.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/15w9ig1258666076.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/16vnt81258666076.tab")
+ }
>
> system("convert tmp/1t1kx1258666076.ps tmp/1t1kx1258666076.png")
> system("convert tmp/2a7zg1258666076.ps tmp/2a7zg1258666076.png")
> system("convert tmp/39onr1258666076.ps tmp/39onr1258666076.png")
> system("convert tmp/49vfk1258666076.ps tmp/49vfk1258666076.png")
> system("convert tmp/5fmmj1258666076.ps tmp/5fmmj1258666076.png")
> system("convert tmp/6kvie1258666076.ps tmp/6kvie1258666076.png")
> system("convert tmp/7w07u1258666076.ps tmp/7w07u1258666076.png")
> system("convert tmp/8wbxj1258666076.ps tmp/8wbxj1258666076.png")
> system("convert tmp/9myqd1258666076.ps tmp/9myqd1258666076.png")
> system("convert tmp/10xqmf1258666076.ps tmp/10xqmf1258666076.png")
>
>
> proc.time()
user system elapsed
2.342 1.602 2.922