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(-1.2
+ ,23.6
+ ,0.2
+ ,-2.2
+ ,-2.4
+ ,25.7
+ ,-1.2
+ ,-4.2
+ ,0.8
+ ,32.5
+ ,-2.4
+ ,-1.6
+ ,-0.1
+ ,33.5
+ ,0.8
+ ,-1.9
+ ,-1.5
+ ,34.5
+ ,-0.1
+ ,0.2
+ ,-4.4
+ ,27.9
+ ,-1.5
+ ,-1.2
+ ,-4.2
+ ,45.3
+ ,-4.4
+ ,-2.4
+ ,3.5
+ ,40.8
+ ,-4.2
+ ,0.8
+ ,10
+ ,58.5
+ ,3.5
+ ,-0.1
+ ,8.6
+ ,32.5
+ ,10
+ ,-1.5
+ ,9.5
+ ,35.5
+ ,8.6
+ ,-4.4
+ ,9.9
+ ,46.7
+ ,9.5
+ ,-4.2
+ ,10.4
+ ,53.2
+ ,9.9
+ ,3.5
+ ,16
+ ,36.1
+ ,10.4
+ ,10
+ ,12.7
+ ,54
+ ,16
+ ,8.6
+ ,10.2
+ ,58.1
+ ,12.7
+ ,9.5
+ ,8.9
+ ,41.8
+ ,10.2
+ ,9.9
+ ,12.6
+ ,43.1
+ ,8.9
+ ,10.4
+ ,13.6
+ ,76
+ ,12.6
+ ,16
+ ,14.8
+ ,42.8
+ ,13.6
+ ,12.7
+ ,9.5
+ ,41
+ ,14.8
+ ,10.2
+ ,13.7
+ ,61.4
+ ,9.5
+ ,8.9
+ ,17
+ ,34.2
+ ,13.7
+ ,12.6
+ ,14.7
+ ,53.8
+ ,17
+ ,13.6
+ ,17.4
+ ,80.7
+ ,14.7
+ ,14.8
+ ,9
+ ,79.5
+ ,17.4
+ ,9.5
+ ,9.1
+ ,96.5
+ ,9
+ ,13.7
+ ,12.2
+ ,108.3
+ ,9.1
+ ,17
+ ,15.9
+ ,100.1
+ ,12.2
+ ,14.7
+ ,12.9
+ ,108.5
+ ,15.9
+ ,17.4
+ ,10.9
+ ,127.4
+ ,12.9
+ ,9
+ ,10.6
+ ,86.5
+ ,10.9
+ ,9.1
+ ,13.2
+ ,71.4
+ ,10.6
+ ,12.2
+ ,9.6
+ ,88.2
+ ,13.2
+ ,15.9
+ ,6.4
+ ,135.6
+ ,9.6
+ ,12.9
+ ,5.8
+ ,70.5
+ ,6.4
+ ,10.9
+ ,-1
+ ,87.5
+ ,5.8
+ ,10.6
+ ,-0.2
+ ,73.3
+ ,-1
+ ,13.2
+ ,2.7
+ ,92.2
+ ,-0.2
+ ,9.6
+ ,3.6
+ ,61.1
+ ,2.7
+ ,6.4
+ ,-0.9
+ ,45.7
+ ,3.6
+ ,5.8
+ ,0.3
+ ,30.5
+ ,-0.9
+ ,-1
+ ,-1.1
+ ,34.8
+ ,0.3
+ ,-0.2
+ ,-2.5
+ ,29.2
+ ,-1.1
+ ,2.7
+ ,-3.4
+ ,56.7
+ ,-2.5
+ ,3.6
+ ,-3.5
+ ,67.1
+ ,-3.4
+ ,-0.9
+ ,-3.9
+ ,41.8
+ ,-3.5
+ ,0.3
+ ,-4.6
+ ,46.8
+ ,-3.9
+ ,-1.1
+ ,-0.1
+ ,50.1
+ ,-4.6
+ ,-2.5
+ ,4.3
+ ,81.9
+ ,-0.1
+ ,-3.4
+ ,10.2
+ ,115.8
+ ,4.3
+ ,-3.5
+ ,8.7
+ ,102.5
+ ,10.2
+ ,-3.9
+ ,13.3
+ ,106.6
+ ,8.7
+ ,-4.6
+ ,15
+ ,101.4
+ ,13.3
+ ,-0.1
+ ,20.7
+ ,136.1
+ ,15
+ ,4.3
+ ,20.7
+ ,143.4
+ ,20.7
+ ,10.2
+ ,26.4
+ ,127.5
+ ,20.7
+ ,8.7
+ ,31.2
+ ,113.8
+ ,26.4
+ ,13.3
+ ,31.4
+ ,75.3
+ ,31.2
+ ,15
+ ,26.6
+ ,98.5
+ ,31.4
+ ,20.7
+ ,26.6
+ ,113.7
+ ,26.6
+ ,20.7
+ ,19.2
+ ,103.7
+ ,26.6
+ ,26.4
+ ,6.5
+ ,73.9
+ ,19.2
+ ,31.2
+ ,3.1
+ ,52.5
+ ,6.5
+ ,31.4
+ ,-0.2
+ ,63.9
+ ,3.1
+ ,26.6
+ ,-4
+ ,44.9
+ ,-0.2
+ ,26.6
+ ,-12.6
+ ,31.3
+ ,-4
+ ,19.2
+ ,-13
+ ,24.9
+ ,-12.6
+ ,6.5
+ ,-17.6
+ ,22.8
+ ,-13
+ ,3.1
+ ,-21.7
+ ,24.8
+ ,-17.6
+ ,-0.2
+ ,-23.2
+ ,22.8
+ ,-21.7
+ ,-4
+ ,-16.8
+ ,20.9
+ ,-23.2
+ ,-12.6
+ ,-19.8
+ ,21.5
+ ,-16.8
+ ,-13)
+ ,dim=c(4
+ ,73)
+ ,dimnames=list(c('Y'
+ ,'X'
+ ,'Y1'
+ ,'Y5')
+ ,1:73))
> y <- array(NA,dim=c(4,73),dimnames=list(c('Y','X','Y1','Y5'),1:73))
> 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
Y X Y1 Y5 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 -1.2 23.6 0.2 -2.2 1 0 0 0 0 0 0 0 0 0 0 1
2 -2.4 25.7 -1.2 -4.2 0 1 0 0 0 0 0 0 0 0 0 2
3 0.8 32.5 -2.4 -1.6 0 0 1 0 0 0 0 0 0 0 0 3
4 -0.1 33.5 0.8 -1.9 0 0 0 1 0 0 0 0 0 0 0 4
5 -1.5 34.5 -0.1 0.2 0 0 0 0 1 0 0 0 0 0 0 5
6 -4.4 27.9 -1.5 -1.2 0 0 0 0 0 1 0 0 0 0 0 6
7 -4.2 45.3 -4.4 -2.4 0 0 0 0 0 0 1 0 0 0 0 7
8 3.5 40.8 -4.2 0.8 0 0 0 0 0 0 0 1 0 0 0 8
9 10.0 58.5 3.5 -0.1 0 0 0 0 0 0 0 0 1 0 0 9
10 8.6 32.5 10.0 -1.5 0 0 0 0 0 0 0 0 0 1 0 10
11 9.5 35.5 8.6 -4.4 0 0 0 0 0 0 0 0 0 0 1 11
12 9.9 46.7 9.5 -4.2 0 0 0 0 0 0 0 0 0 0 0 12
13 10.4 53.2 9.9 3.5 1 0 0 0 0 0 0 0 0 0 0 13
14 16.0 36.1 10.4 10.0 0 1 0 0 0 0 0 0 0 0 0 14
15 12.7 54.0 16.0 8.6 0 0 1 0 0 0 0 0 0 0 0 15
16 10.2 58.1 12.7 9.5 0 0 0 1 0 0 0 0 0 0 0 16
17 8.9 41.8 10.2 9.9 0 0 0 0 1 0 0 0 0 0 0 17
18 12.6 43.1 8.9 10.4 0 0 0 0 0 1 0 0 0 0 0 18
19 13.6 76.0 12.6 16.0 0 0 0 0 0 0 1 0 0 0 0 19
20 14.8 42.8 13.6 12.7 0 0 0 0 0 0 0 1 0 0 0 20
21 9.5 41.0 14.8 10.2 0 0 0 0 0 0 0 0 1 0 0 21
22 13.7 61.4 9.5 8.9 0 0 0 0 0 0 0 0 0 1 0 22
23 17.0 34.2 13.7 12.6 0 0 0 0 0 0 0 0 0 0 1 23
24 14.7 53.8 17.0 13.6 0 0 0 0 0 0 0 0 0 0 0 24
25 17.4 80.7 14.7 14.8 1 0 0 0 0 0 0 0 0 0 0 25
26 9.0 79.5 17.4 9.5 0 1 0 0 0 0 0 0 0 0 0 26
27 9.1 96.5 9.0 13.7 0 0 1 0 0 0 0 0 0 0 0 27
28 12.2 108.3 9.1 17.0 0 0 0 1 0 0 0 0 0 0 0 28
29 15.9 100.1 12.2 14.7 0 0 0 0 1 0 0 0 0 0 0 29
30 12.9 108.5 15.9 17.4 0 0 0 0 0 1 0 0 0 0 0 30
31 10.9 127.4 12.9 9.0 0 0 0 0 0 0 1 0 0 0 0 31
32 10.6 86.5 10.9 9.1 0 0 0 0 0 0 0 1 0 0 0 32
33 13.2 71.4 10.6 12.2 0 0 0 0 0 0 0 0 1 0 0 33
34 9.6 88.2 13.2 15.9 0 0 0 0 0 0 0 0 0 1 0 34
35 6.4 135.6 9.6 12.9 0 0 0 0 0 0 0 0 0 0 1 35
36 5.8 70.5 6.4 10.9 0 0 0 0 0 0 0 0 0 0 0 36
37 -1.0 87.5 5.8 10.6 1 0 0 0 0 0 0 0 0 0 0 37
38 -0.2 73.3 -1.0 13.2 0 1 0 0 0 0 0 0 0 0 0 38
39 2.7 92.2 -0.2 9.6 0 0 1 0 0 0 0 0 0 0 0 39
40 3.6 61.1 2.7 6.4 0 0 0 1 0 0 0 0 0 0 0 40
41 -0.9 45.7 3.6 5.8 0 0 0 0 1 0 0 0 0 0 0 41
42 0.3 30.5 -0.9 -1.0 0 0 0 0 0 1 0 0 0 0 0 42
43 -1.1 34.8 0.3 -0.2 0 0 0 0 0 0 1 0 0 0 0 43
44 -2.5 29.2 -1.1 2.7 0 0 0 0 0 0 0 1 0 0 0 44
45 -3.4 56.7 -2.5 3.6 0 0 0 0 0 0 0 0 1 0 0 45
46 -3.5 67.1 -3.4 -0.9 0 0 0 0 0 0 0 0 0 1 0 46
47 -3.9 41.8 -3.5 0.3 0 0 0 0 0 0 0 0 0 0 1 47
48 -4.6 46.8 -3.9 -1.1 0 0 0 0 0 0 0 0 0 0 0 48
49 -0.1 50.1 -4.6 -2.5 1 0 0 0 0 0 0 0 0 0 0 49
50 4.3 81.9 -0.1 -3.4 0 1 0 0 0 0 0 0 0 0 0 50
51 10.2 115.8 4.3 -3.5 0 0 1 0 0 0 0 0 0 0 0 51
52 8.7 102.5 10.2 -3.9 0 0 0 1 0 0 0 0 0 0 0 52
53 13.3 106.6 8.7 -4.6 0 0 0 0 1 0 0 0 0 0 0 53
54 15.0 101.4 13.3 -0.1 0 0 0 0 0 1 0 0 0 0 0 54
55 20.7 136.1 15.0 4.3 0 0 0 0 0 0 1 0 0 0 0 55
56 20.7 143.4 20.7 10.2 0 0 0 0 0 0 0 1 0 0 0 56
57 26.4 127.5 20.7 8.7 0 0 0 0 0 0 0 0 1 0 0 57
58 31.2 113.8 26.4 13.3 0 0 0 0 0 0 0 0 0 1 0 58
59 31.4 75.3 31.2 15.0 0 0 0 0 0 0 0 0 0 0 1 59
60 26.6 98.5 31.4 20.7 0 0 0 0 0 0 0 0 0 0 0 60
61 26.6 113.7 26.6 20.7 1 0 0 0 0 0 0 0 0 0 0 61
62 19.2 103.7 26.6 26.4 0 1 0 0 0 0 0 0 0 0 0 62
63 6.5 73.9 19.2 31.2 0 0 1 0 0 0 0 0 0 0 0 63
64 3.1 52.5 6.5 31.4 0 0 0 1 0 0 0 0 0 0 0 64
65 -0.2 63.9 3.1 26.6 0 0 0 0 1 0 0 0 0 0 0 65
66 -4.0 44.9 -0.2 26.6 0 0 0 0 0 1 0 0 0 0 0 66
67 -12.6 31.3 -4.0 19.2 0 0 0 0 0 0 1 0 0 0 0 67
68 -13.0 24.9 -12.6 6.5 0 0 0 0 0 0 0 1 0 0 0 68
69 -17.6 22.8 -13.0 3.1 0 0 0 0 0 0 0 0 1 0 0 69
70 -21.7 24.8 -17.6 -0.2 0 0 0 0 0 0 0 0 0 1 0 70
71 -23.2 22.8 -21.7 -4.0 0 0 0 0 0 0 0 0 0 0 1 71
72 -16.8 20.9 -23.2 -12.6 0 0 0 0 0 0 0 0 0 0 0 72
73 -19.8 21.5 -16.8 -13.0 1 0 0 0 0 0 0 0 0 0 0 73
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X Y1 Y5 M1 M2
-0.228904 0.051222 0.953264 -0.166874 -0.775207 -0.959995
M3 M4 M5 M6 M7 M8
-0.952058 -0.573275 -0.176635 -0.003322 -1.296379 1.294355
M9 M10 M11 t
0.717593 -0.052491 0.184538 -0.044462
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-7.7287 -2.3454 -0.2636 2.7361 7.7171
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.228904 1.840167 -0.124 0.9014
X 0.051222 0.019664 2.605 0.0117 *
Y1 0.953264 0.068625 13.891 <2e-16 ***
Y5 -0.166874 0.062741 -2.660 0.0101 *
M1 -0.775207 2.075813 -0.373 0.7102
M2 -0.959995 2.180140 -0.440 0.6614
M3 -0.952058 2.237315 -0.426 0.6720
M4 -0.573275 2.208944 -0.260 0.7962
M5 -0.176635 2.189287 -0.081 0.9360
M6 -0.003322 2.173918 -0.002 0.9988
M7 -1.296379 2.216391 -0.585 0.5609
M8 1.294355 2.169092 0.597 0.5531
M9 0.717593 2.154867 0.333 0.7403
M10 -0.052491 2.149686 -0.024 0.9806
M11 0.184538 2.138932 0.086 0.9315
t -0.044462 0.025632 -1.735 0.0882 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.702 on 57 degrees of freedom
Multiple R-squared: 0.9207, Adjusted R-squared: 0.8998
F-statistic: 44.13 on 15 and 57 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.2303219 0.4606437 0.7696781
[2,] 0.1152133 0.2304265 0.8847867
[3,] 0.3230686 0.6461373 0.6769314
[4,] 0.4157363 0.8314727 0.5842637
[5,] 0.3438586 0.6877172 0.6561414
[6,] 0.3379177 0.6758355 0.6620823
[7,] 0.2938811 0.5877622 0.7061189
[8,] 0.5572965 0.8854071 0.4427035
[9,] 0.4821026 0.9642052 0.5178974
[10,] 0.4055458 0.8110915 0.5944542
[11,] 0.3986154 0.7972309 0.6013846
[12,] 0.4282884 0.8565768 0.5717116
[13,] 0.3662249 0.7324498 0.6337751
[14,] 0.2820556 0.5641111 0.7179444
[15,] 0.2580438 0.5160876 0.7419562
[16,] 0.3145028 0.6290057 0.6854972
[17,] 0.5053155 0.9893691 0.4946845
[18,] 0.4183635 0.8367270 0.5816365
[19,] 0.5791230 0.8417539 0.4208770
[20,] 0.5026538 0.9946924 0.4973462
[21,] 0.4840681 0.9681362 0.5159319
[22,] 0.4436473 0.8872945 0.5563527
[23,] 0.3918885 0.7837771 0.6081115
[24,] 0.3639357 0.7278715 0.6360643
[25,] 0.2849856 0.5699713 0.7150144
[26,] 0.2342984 0.4685968 0.7657016
[27,] 0.1760935 0.3521871 0.8239065
[28,] 0.1438111 0.2876221 0.8561889
[29,] 0.1108452 0.2216905 0.8891548
[30,] 0.1316244 0.2632487 0.8683756
[31,] 0.1371126 0.2742252 0.8628874
[32,] 0.1212613 0.2425226 0.8787387
[33,] 0.1536606 0.3073213 0.8463394
[34,] 0.1431248 0.2862496 0.8568752
[35,] 0.0926310 0.1852620 0.9073690
[36,] 0.1682988 0.3365975 0.8317012
> postscript(file="/var/www/html/rcomp/tmp/1biib1261321821.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/2pmcu1261321821.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/3rn1w1261321821.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/44zlp1261321821.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/5lofz1261321821.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 = 73
Frequency = 1
1 2 3 4 5 6 7
-1.9180451 -1.9955413 2.4704634 -1.9155851 -2.5106121 -4.1004523 -0.8899825
8 9 10 11 12 13 14
4.8375902 3.5618664 -2.1216512 -0.7172513 -1.4865010 0.4038518 7.7170512
15 16 17 18 19 20 21
-2.0351990 -1.7835734 -0.1509234 4.6763171 2.7360514 1.5864028 -4.5612766
22 23 24 25 26 27 28
4.2436995 5.3581012 -0.6957470 3.8388034 -7.7287263 0.2453102 2.8619283
29 30 31 32 33 34 35
3.2908426 -3.3447876 -3.5173206 -2.3453954 2.4525709 -3.0544636 -5.9438311
36 37 38 39 40 41 42
-0.2635802 -6.5927910 2.0798782 2.6849472 1.5451702 -3.4762502 1.5284145
43 44 45 46 47 48 49
0.2352616 -1.6056620 -1.8082896 -1.5194508 -0.5205238 -1.0999532 4.4843428
50 51 52 53 54 55 56
3.0448579 3.0339076 -3.8101644 1.5407299 -0.2558434 4.1179690 -3.2512673
57 58 59 60 61 62 63
3.6340752 5.2843837 2.9718877 -2.0269329 2.5898253 -3.1175198 -6.3994294
64 65 66 67 68 69 70
3.1022244 1.3062132 1.4963517 -2.6819788 0.7783318 -3.2789463 -2.8325176
71 72 73
-1.1483827 5.5727142 -2.8059872
> postscript(file="/var/www/html/rcomp/tmp/6v8x31261321821.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 = 73
Frequency = 1
lag(myerror, k = 1) myerror
0 -1.9180451 NA
1 -1.9955413 -1.9180451
2 2.4704634 -1.9955413
3 -1.9155851 2.4704634
4 -2.5106121 -1.9155851
5 -4.1004523 -2.5106121
6 -0.8899825 -4.1004523
7 4.8375902 -0.8899825
8 3.5618664 4.8375902
9 -2.1216512 3.5618664
10 -0.7172513 -2.1216512
11 -1.4865010 -0.7172513
12 0.4038518 -1.4865010
13 7.7170512 0.4038518
14 -2.0351990 7.7170512
15 -1.7835734 -2.0351990
16 -0.1509234 -1.7835734
17 4.6763171 -0.1509234
18 2.7360514 4.6763171
19 1.5864028 2.7360514
20 -4.5612766 1.5864028
21 4.2436995 -4.5612766
22 5.3581012 4.2436995
23 -0.6957470 5.3581012
24 3.8388034 -0.6957470
25 -7.7287263 3.8388034
26 0.2453102 -7.7287263
27 2.8619283 0.2453102
28 3.2908426 2.8619283
29 -3.3447876 3.2908426
30 -3.5173206 -3.3447876
31 -2.3453954 -3.5173206
32 2.4525709 -2.3453954
33 -3.0544636 2.4525709
34 -5.9438311 -3.0544636
35 -0.2635802 -5.9438311
36 -6.5927910 -0.2635802
37 2.0798782 -6.5927910
38 2.6849472 2.0798782
39 1.5451702 2.6849472
40 -3.4762502 1.5451702
41 1.5284145 -3.4762502
42 0.2352616 1.5284145
43 -1.6056620 0.2352616
44 -1.8082896 -1.6056620
45 -1.5194508 -1.8082896
46 -0.5205238 -1.5194508
47 -1.0999532 -0.5205238
48 4.4843428 -1.0999532
49 3.0448579 4.4843428
50 3.0339076 3.0448579
51 -3.8101644 3.0339076
52 1.5407299 -3.8101644
53 -0.2558434 1.5407299
54 4.1179690 -0.2558434
55 -3.2512673 4.1179690
56 3.6340752 -3.2512673
57 5.2843837 3.6340752
58 2.9718877 5.2843837
59 -2.0269329 2.9718877
60 2.5898253 -2.0269329
61 -3.1175198 2.5898253
62 -6.3994294 -3.1175198
63 3.1022244 -6.3994294
64 1.3062132 3.1022244
65 1.4963517 1.3062132
66 -2.6819788 1.4963517
67 0.7783318 -2.6819788
68 -3.2789463 0.7783318
69 -2.8325176 -3.2789463
70 -1.1483827 -2.8325176
71 5.5727142 -1.1483827
72 -2.8059872 5.5727142
73 NA -2.8059872
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1.9955413 -1.9180451
[2,] 2.4704634 -1.9955413
[3,] -1.9155851 2.4704634
[4,] -2.5106121 -1.9155851
[5,] -4.1004523 -2.5106121
[6,] -0.8899825 -4.1004523
[7,] 4.8375902 -0.8899825
[8,] 3.5618664 4.8375902
[9,] -2.1216512 3.5618664
[10,] -0.7172513 -2.1216512
[11,] -1.4865010 -0.7172513
[12,] 0.4038518 -1.4865010
[13,] 7.7170512 0.4038518
[14,] -2.0351990 7.7170512
[15,] -1.7835734 -2.0351990
[16,] -0.1509234 -1.7835734
[17,] 4.6763171 -0.1509234
[18,] 2.7360514 4.6763171
[19,] 1.5864028 2.7360514
[20,] -4.5612766 1.5864028
[21,] 4.2436995 -4.5612766
[22,] 5.3581012 4.2436995
[23,] -0.6957470 5.3581012
[24,] 3.8388034 -0.6957470
[25,] -7.7287263 3.8388034
[26,] 0.2453102 -7.7287263
[27,] 2.8619283 0.2453102
[28,] 3.2908426 2.8619283
[29,] -3.3447876 3.2908426
[30,] -3.5173206 -3.3447876
[31,] -2.3453954 -3.5173206
[32,] 2.4525709 -2.3453954
[33,] -3.0544636 2.4525709
[34,] -5.9438311 -3.0544636
[35,] -0.2635802 -5.9438311
[36,] -6.5927910 -0.2635802
[37,] 2.0798782 -6.5927910
[38,] 2.6849472 2.0798782
[39,] 1.5451702 2.6849472
[40,] -3.4762502 1.5451702
[41,] 1.5284145 -3.4762502
[42,] 0.2352616 1.5284145
[43,] -1.6056620 0.2352616
[44,] -1.8082896 -1.6056620
[45,] -1.5194508 -1.8082896
[46,] -0.5205238 -1.5194508
[47,] -1.0999532 -0.5205238
[48,] 4.4843428 -1.0999532
[49,] 3.0448579 4.4843428
[50,] 3.0339076 3.0448579
[51,] -3.8101644 3.0339076
[52,] 1.5407299 -3.8101644
[53,] -0.2558434 1.5407299
[54,] 4.1179690 -0.2558434
[55,] -3.2512673 4.1179690
[56,] 3.6340752 -3.2512673
[57,] 5.2843837 3.6340752
[58,] 2.9718877 5.2843837
[59,] -2.0269329 2.9718877
[60,] 2.5898253 -2.0269329
[61,] -3.1175198 2.5898253
[62,] -6.3994294 -3.1175198
[63,] 3.1022244 -6.3994294
[64,] 1.3062132 3.1022244
[65,] 1.4963517 1.3062132
[66,] -2.6819788 1.4963517
[67,] 0.7783318 -2.6819788
[68,] -3.2789463 0.7783318
[69,] -2.8325176 -3.2789463
[70,] -1.1483827 -2.8325176
[71,] 5.5727142 -1.1483827
[72,] -2.8059872 5.5727142
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1.9955413 -1.9180451
2 2.4704634 -1.9955413
3 -1.9155851 2.4704634
4 -2.5106121 -1.9155851
5 -4.1004523 -2.5106121
6 -0.8899825 -4.1004523
7 4.8375902 -0.8899825
8 3.5618664 4.8375902
9 -2.1216512 3.5618664
10 -0.7172513 -2.1216512
11 -1.4865010 -0.7172513
12 0.4038518 -1.4865010
13 7.7170512 0.4038518
14 -2.0351990 7.7170512
15 -1.7835734 -2.0351990
16 -0.1509234 -1.7835734
17 4.6763171 -0.1509234
18 2.7360514 4.6763171
19 1.5864028 2.7360514
20 -4.5612766 1.5864028
21 4.2436995 -4.5612766
22 5.3581012 4.2436995
23 -0.6957470 5.3581012
24 3.8388034 -0.6957470
25 -7.7287263 3.8388034
26 0.2453102 -7.7287263
27 2.8619283 0.2453102
28 3.2908426 2.8619283
29 -3.3447876 3.2908426
30 -3.5173206 -3.3447876
31 -2.3453954 -3.5173206
32 2.4525709 -2.3453954
33 -3.0544636 2.4525709
34 -5.9438311 -3.0544636
35 -0.2635802 -5.9438311
36 -6.5927910 -0.2635802
37 2.0798782 -6.5927910
38 2.6849472 2.0798782
39 1.5451702 2.6849472
40 -3.4762502 1.5451702
41 1.5284145 -3.4762502
42 0.2352616 1.5284145
43 -1.6056620 0.2352616
44 -1.8082896 -1.6056620
45 -1.5194508 -1.8082896
46 -0.5205238 -1.5194508
47 -1.0999532 -0.5205238
48 4.4843428 -1.0999532
49 3.0448579 4.4843428
50 3.0339076 3.0448579
51 -3.8101644 3.0339076
52 1.5407299 -3.8101644
53 -0.2558434 1.5407299
54 4.1179690 -0.2558434
55 -3.2512673 4.1179690
56 3.6340752 -3.2512673
57 5.2843837 3.6340752
58 2.9718877 5.2843837
59 -2.0269329 2.9718877
60 2.5898253 -2.0269329
61 -3.1175198 2.5898253
62 -6.3994294 -3.1175198
63 3.1022244 -6.3994294
64 1.3062132 3.1022244
65 1.4963517 1.3062132
66 -2.6819788 1.4963517
67 0.7783318 -2.6819788
68 -3.2789463 0.7783318
69 -2.8325176 -3.2789463
70 -1.1483827 -2.8325176
71 5.5727142 -1.1483827
72 -2.8059872 5.5727142
> 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/7eeey1261321821.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/8u89d1261321821.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/9lzmi1261321821.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/10lzoa1261321821.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/110se71261321821.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/129qbs1261321821.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/13ysbw1261321821.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/149j661261321822.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/1519i81261321822.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/16ziqd1261321822.tab")
+ }
>
> try(system("convert tmp/1biib1261321821.ps tmp/1biib1261321821.png",intern=TRUE))
character(0)
> try(system("convert tmp/2pmcu1261321821.ps tmp/2pmcu1261321821.png",intern=TRUE))
character(0)
> try(system("convert tmp/3rn1w1261321821.ps tmp/3rn1w1261321821.png",intern=TRUE))
character(0)
> try(system("convert tmp/44zlp1261321821.ps tmp/44zlp1261321821.png",intern=TRUE))
character(0)
> try(system("convert tmp/5lofz1261321821.ps tmp/5lofz1261321821.png",intern=TRUE))
character(0)
> try(system("convert tmp/6v8x31261321821.ps tmp/6v8x31261321821.png",intern=TRUE))
character(0)
> try(system("convert tmp/7eeey1261321821.ps tmp/7eeey1261321821.png",intern=TRUE))
character(0)
> try(system("convert tmp/8u89d1261321821.ps tmp/8u89d1261321821.png",intern=TRUE))
character(0)
> try(system("convert tmp/9lzmi1261321821.ps tmp/9lzmi1261321821.png",intern=TRUE))
character(0)
> try(system("convert tmp/10lzoa1261321821.ps tmp/10lzoa1261321821.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.562 1.579 3.399