R version 2.12.0 (2010-10-15)
Copyright (C) 2010 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-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(29462
+ ,27071
+ ,31514
+ ,26105
+ ,29462
+ ,27071
+ ,22397
+ ,26105
+ ,29462
+ ,23843
+ ,22397
+ ,26105
+ ,21705
+ ,23843
+ ,22397
+ ,18089
+ ,21705
+ ,23843
+ ,20764
+ ,18089
+ ,21705
+ ,25316
+ ,20764
+ ,18089
+ ,17704
+ ,25316
+ ,20764
+ ,15548
+ ,17704
+ ,25316
+ ,28029
+ ,15548
+ ,17704
+ ,29383
+ ,28029
+ ,15548
+ ,36438
+ ,29383
+ ,28029
+ ,32034
+ ,36438
+ ,29383
+ ,22679
+ ,32034
+ ,36438
+ ,24319
+ ,22679
+ ,32034
+ ,18004
+ ,24319
+ ,22679
+ ,17537
+ ,18004
+ ,24319
+ ,20366
+ ,17537
+ ,18004
+ ,22782
+ ,20366
+ ,17537
+ ,19169
+ ,22782
+ ,20366
+ ,13807
+ ,19169
+ ,22782
+ ,29743
+ ,13807
+ ,19169
+ ,25591
+ ,29743
+ ,13807
+ ,29096
+ ,25591
+ ,29743
+ ,26482
+ ,29096
+ ,25591
+ ,22405
+ ,26482
+ ,29096
+ ,27044
+ ,22405
+ ,26482
+ ,17970
+ ,27044
+ ,22405
+ ,18730
+ ,17970
+ ,27044
+ ,19684
+ ,18730
+ ,17970
+ ,19785
+ ,19684
+ ,18730
+ ,18479
+ ,19785
+ ,19684
+ ,10698
+ ,18479
+ ,19785
+ ,31956
+ ,10698
+ ,18479
+ ,29506
+ ,31956
+ ,10698
+ ,34506
+ ,29506
+ ,31956
+ ,27165
+ ,34506
+ ,29506
+ ,26736
+ ,27165
+ ,34506
+ ,23691
+ ,26736
+ ,27165
+ ,18157
+ ,23691
+ ,26736
+ ,17328
+ ,18157
+ ,23691
+ ,18205
+ ,17328
+ ,18157
+ ,20995
+ ,18205
+ ,17328
+ ,17382
+ ,20995
+ ,18205
+ ,9367
+ ,17382
+ ,20995
+ ,31124
+ ,9367
+ ,17382
+ ,26551
+ ,31124
+ ,9367
+ ,30651
+ ,26551
+ ,31124
+ ,25859
+ ,30651
+ ,26551
+ ,25100
+ ,25859
+ ,30651
+ ,25778
+ ,25100
+ ,25859
+ ,20418
+ ,25778
+ ,25100
+ ,18688
+ ,20418
+ ,25778
+ ,20424
+ ,18688
+ ,20418
+ ,24776
+ ,20424
+ ,18688
+ ,19814
+ ,24776
+ ,20424
+ ,12738
+ ,19814
+ ,24776
+ ,31566
+ ,12738
+ ,19814
+ ,30111
+ ,31566
+ ,12738
+ ,30019
+ ,30111
+ ,31566
+ ,31934
+ ,30019
+ ,30111
+ ,25826
+ ,31934
+ ,30019
+ ,26835
+ ,25826
+ ,31934
+ ,20205
+ ,26835
+ ,25826
+ ,17789
+ ,20205
+ ,26835
+ ,20520
+ ,17789
+ ,20205
+ ,22518
+ ,20520
+ ,17789
+ ,15572
+ ,22518
+ ,20520
+ ,11509
+ ,15572
+ ,22518
+ ,25447
+ ,11509
+ ,15572
+ ,24090
+ ,25447
+ ,11509
+ ,27786
+ ,24090
+ ,25447
+ ,26195
+ ,27786
+ ,24090
+ ,20516
+ ,26195
+ ,27786
+ ,22759
+ ,20516
+ ,26195
+ ,19028
+ ,22759
+ ,20516
+ ,16971
+ ,19028
+ ,22759
+ ,20036
+ ,16971
+ ,19028
+ ,22485
+ ,20036
+ ,16971
+ ,18730
+ ,22485
+ ,20036
+ ,14538
+ ,18730
+ ,22485
+ ,27561
+ ,14538
+ ,18730
+ ,25985
+ ,27561
+ ,14538
+ ,34670
+ ,25985
+ ,27561
+ ,32066
+ ,34670
+ ,25985
+ ,27186
+ ,32066
+ ,34670
+ ,29586
+ ,27186
+ ,32066
+ ,21359
+ ,29586
+ ,27186
+ ,21553
+ ,21359
+ ,29586
+ ,19573
+ ,21553
+ ,21359
+ ,24256
+ ,19573
+ ,21553)
+ ,dim=c(3
+ ,92)
+ ,dimnames=list(c('X'
+ ,'Y_1'
+ ,'Y_2')
+ ,1:92))
> y <- array(NA,dim=c(3,92),dimnames=list(c('X','Y_1','Y_2'),1:92))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
> 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
X Y_1 Y_2 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 29462 27071 31514 1 0 0 0 0 0 0 0 0 0 0
2 26105 29462 27071 0 1 0 0 0 0 0 0 0 0 0
3 22397 26105 29462 0 0 1 0 0 0 0 0 0 0 0
4 23843 22397 26105 0 0 0 1 0 0 0 0 0 0 0
5 21705 23843 22397 0 0 0 0 1 0 0 0 0 0 0
6 18089 21705 23843 0 0 0 0 0 1 0 0 0 0 0
7 20764 18089 21705 0 0 0 0 0 0 1 0 0 0 0
8 25316 20764 18089 0 0 0 0 0 0 0 1 0 0 0
9 17704 25316 20764 0 0 0 0 0 0 0 0 1 0 0
10 15548 17704 25316 0 0 0 0 0 0 0 0 0 1 0
11 28029 15548 17704 0 0 0 0 0 0 0 0 0 0 1
12 29383 28029 15548 0 0 0 0 0 0 0 0 0 0 0
13 36438 29383 28029 1 0 0 0 0 0 0 0 0 0 0
14 32034 36438 29383 0 1 0 0 0 0 0 0 0 0 0
15 22679 32034 36438 0 0 1 0 0 0 0 0 0 0 0
16 24319 22679 32034 0 0 0 1 0 0 0 0 0 0 0
17 18004 24319 22679 0 0 0 0 1 0 0 0 0 0 0
18 17537 18004 24319 0 0 0 0 0 1 0 0 0 0 0
19 20366 17537 18004 0 0 0 0 0 0 1 0 0 0 0
20 22782 20366 17537 0 0 0 0 0 0 0 1 0 0 0
21 19169 22782 20366 0 0 0 0 0 0 0 0 1 0 0
22 13807 19169 22782 0 0 0 0 0 0 0 0 0 1 0
23 29743 13807 19169 0 0 0 0 0 0 0 0 0 0 1
24 25591 29743 13807 0 0 0 0 0 0 0 0 0 0 0
25 29096 25591 29743 1 0 0 0 0 0 0 0 0 0 0
26 26482 29096 25591 0 1 0 0 0 0 0 0 0 0 0
27 22405 26482 29096 0 0 1 0 0 0 0 0 0 0 0
28 27044 22405 26482 0 0 0 1 0 0 0 0 0 0 0
29 17970 27044 22405 0 0 0 0 1 0 0 0 0 0 0
30 18730 17970 27044 0 0 0 0 0 1 0 0 0 0 0
31 19684 18730 17970 0 0 0 0 0 0 1 0 0 0 0
32 19785 19684 18730 0 0 0 0 0 0 0 1 0 0 0
33 18479 19785 19684 0 0 0 0 0 0 0 0 1 0 0
34 10698 18479 19785 0 0 0 0 0 0 0 0 0 1 0
35 31956 10698 18479 0 0 0 0 0 0 0 0 0 0 1
36 29506 31956 10698 0 0 0 0 0 0 0 0 0 0 0
37 34506 29506 31956 1 0 0 0 0 0 0 0 0 0 0
38 27165 34506 29506 0 1 0 0 0 0 0 0 0 0 0
39 26736 27165 34506 0 0 1 0 0 0 0 0 0 0 0
40 23691 26736 27165 0 0 0 1 0 0 0 0 0 0 0
41 18157 23691 26736 0 0 0 0 1 0 0 0 0 0 0
42 17328 18157 23691 0 0 0 0 0 1 0 0 0 0 0
43 18205 17328 18157 0 0 0 0 0 0 1 0 0 0 0
44 20995 18205 17328 0 0 0 0 0 0 0 1 0 0 0
45 17382 20995 18205 0 0 0 0 0 0 0 0 1 0 0
46 9367 17382 20995 0 0 0 0 0 0 0 0 0 1 0
47 31124 9367 17382 0 0 0 0 0 0 0 0 0 0 1
48 26551 31124 9367 0 0 0 0 0 0 0 0 0 0 0
49 30651 26551 31124 1 0 0 0 0 0 0 0 0 0 0
50 25859 30651 26551 0 1 0 0 0 0 0 0 0 0 0
51 25100 25859 30651 0 0 1 0 0 0 0 0 0 0 0
52 25778 25100 25859 0 0 0 1 0 0 0 0 0 0 0
53 20418 25778 25100 0 0 0 0 1 0 0 0 0 0 0
54 18688 20418 25778 0 0 0 0 0 1 0 0 0 0 0
55 20424 18688 20418 0 0 0 0 0 0 1 0 0 0 0
56 24776 20424 18688 0 0 0 0 0 0 0 1 0 0 0
57 19814 24776 20424 0 0 0 0 0 0 0 0 1 0 0
58 12738 19814 24776 0 0 0 0 0 0 0 0 0 1 0
59 31566 12738 19814 0 0 0 0 0 0 0 0 0 0 1
60 30111 31566 12738 0 0 0 0 0 0 0 0 0 0 0
61 30019 30111 31566 1 0 0 0 0 0 0 0 0 0 0
62 31934 30019 30111 0 1 0 0 0 0 0 0 0 0 0
63 25826 31934 30019 0 0 1 0 0 0 0 0 0 0 0
64 26835 25826 31934 0 0 0 1 0 0 0 0 0 0 0
65 20205 26835 25826 0 0 0 0 1 0 0 0 0 0 0
66 17789 20205 26835 0 0 0 0 0 1 0 0 0 0 0
67 20520 17789 20205 0 0 0 0 0 0 1 0 0 0 0
68 22518 20520 17789 0 0 0 0 0 0 0 1 0 0 0
69 15572 22518 20520 0 0 0 0 0 0 0 0 1 0 0
70 11509 15572 22518 0 0 0 0 0 0 0 0 0 1 0
71 25447 11509 15572 0 0 0 0 0 0 0 0 0 0 1
72 24090 25447 11509 0 0 0 0 0 0 0 0 0 0 0
73 27786 24090 25447 1 0 0 0 0 0 0 0 0 0 0
74 26195 27786 24090 0 1 0 0 0 0 0 0 0 0 0
75 20516 26195 27786 0 0 1 0 0 0 0 0 0 0 0
76 22759 20516 26195 0 0 0 1 0 0 0 0 0 0 0
77 19028 22759 20516 0 0 0 0 1 0 0 0 0 0 0
78 16971 19028 22759 0 0 0 0 0 1 0 0 0 0 0
79 20036 16971 19028 0 0 0 0 0 0 1 0 0 0 0
80 22485 20036 16971 0 0 0 0 0 0 0 1 0 0 0
81 18730 22485 20036 0 0 0 0 0 0 0 0 1 0 0
82 14538 18730 22485 0 0 0 0 0 0 0 0 0 1 0
83 27561 14538 18730 0 0 0 0 0 0 0 0 0 0 1
84 25985 27561 14538 0 0 0 0 0 0 0 0 0 0 0
85 34670 25985 27561 1 0 0 0 0 0 0 0 0 0 0
86 32066 34670 25985 0 1 0 0 0 0 0 0 0 0 0
87 27186 32066 34670 0 0 1 0 0 0 0 0 0 0 0
88 29586 27186 32066 0 0 0 1 0 0 0 0 0 0 0
89 21359 29586 27186 0 0 0 0 1 0 0 0 0 0 0
90 21553 21359 29586 0 0 0 0 0 1 0 0 0 0 0
91 19573 21553 21359 0 0 0 0 0 0 1 0 0 0 0
92 24256 19573 21553 0 0 0 0 0 0 0 1 0 0 0
> 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
1.535e+04 2.891e-01 2.766e-01 1.499e+02 -3.545e+03 -8.211e+03
M4 M5 M6 M7 M8 M9
-4.713e+03 -9.776e+03 -9.729e+03 -6.125e+03 -3.322e+03 -9.311e+03
M10 M11
-1.426e+04 5.344e+03
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-3798.4 -1362.4 -96.6 1221.3 4692.7
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.535e+04 3.028e+03 5.069 2.63e-06 ***
Y_1 2.891e-01 1.091e-01 2.651 0.009725 **
Y_2 2.766e-01 1.085e-01 2.550 0.012738 *
M1 1.499e+02 2.189e+03 0.068 0.945592
M2 -3.545e+03 1.815e+03 -1.953 0.054390 .
M3 -8.211e+03 2.325e+03 -3.531 0.000699 ***
M4 -4.713e+03 2.255e+03 -2.090 0.039858 *
M5 -9.776e+03 1.777e+03 -5.501 4.64e-07 ***
M6 -9.729e+03 2.295e+03 -4.239 6.12e-05 ***
M7 -6.125e+03 1.937e+03 -3.163 0.002229 **
M8 -3.322e+03 1.718e+03 -1.934 0.056773 .
M9 -9.311e+03 1.647e+03 -5.653 2.48e-07 ***
M10 -1.426e+04 2.197e+03 -6.488 7.28e-09 ***
M11 5.344e+03 2.374e+03 2.251 0.027217 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1936 on 78 degrees of freedom
Multiple R-squared: 0.9004, Adjusted R-squared: 0.8839
F-statistic: 54.27 on 13 and 78 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.8727796 0.2544408 0.1272204
[2,] 0.8546674 0.2906652 0.1453326
[3,] 0.7903198 0.4193604 0.2096802
[4,] 0.7539953 0.4920094 0.2460047
[5,] 0.7586262 0.4827477 0.2413738
[6,] 0.7778517 0.4442966 0.2221483
[7,] 0.7682786 0.4634427 0.2317214
[8,] 0.8548243 0.2903513 0.1451757
[9,] 0.8546659 0.2906682 0.1453341
[10,] 0.8045565 0.3908870 0.1954435
[11,] 0.7481330 0.5037339 0.2518670
[12,] 0.7643479 0.4713043 0.2356521
[13,] 0.7862903 0.4274193 0.2137097
[14,] 0.7510181 0.4979639 0.2489819
[15,] 0.6980588 0.6038825 0.3019412
[16,] 0.7805907 0.4388187 0.2194093
[17,] 0.7440283 0.5119433 0.2559717
[18,] 0.7799354 0.4401291 0.2200646
[19,] 0.8426343 0.3147314 0.1573657
[20,] 0.8288017 0.3423967 0.1711983
[21,] 0.8151930 0.3696141 0.1848070
[22,] 0.8468702 0.3062597 0.1531298
[23,] 0.8767218 0.2465564 0.1232782
[24,] 0.8739409 0.2521183 0.1260591
[25,] 0.8572764 0.2854472 0.1427236
[26,] 0.8138963 0.3722074 0.1861037
[27,] 0.7796407 0.4407186 0.2203593
[28,] 0.7409802 0.5180396 0.2590198
[29,] 0.6854690 0.6290620 0.3145310
[30,] 0.7241069 0.5517862 0.2758931
[31,] 0.7988972 0.4022056 0.2011028
[32,] 0.7500357 0.4999286 0.2499643
[33,] 0.7117099 0.5765801 0.2882901
[34,] 0.7768499 0.4463001 0.2231501
[35,] 0.7803739 0.4392522 0.2196261
[36,] 0.7299933 0.5400134 0.2700067
[37,] 0.6717980 0.6564040 0.3282020
[38,] 0.6029745 0.7940510 0.3970255
[39,] 0.5312214 0.9375572 0.4687786
[40,] 0.5059951 0.9880098 0.4940049
[41,] 0.4587846 0.9175692 0.5412154
[42,] 0.4236197 0.8472395 0.5763803
[43,] 0.5279783 0.9440434 0.4720217
[44,] 0.5595574 0.8808853 0.4404426
[45,] 0.7991859 0.4016281 0.2008141
[46,] 0.8129533 0.3740934 0.1870467
[47,] 0.7689291 0.4621418 0.2310709
[48,] 0.7185816 0.5628369 0.2814184
[49,] 0.6529594 0.6940812 0.3470406
[50,] 0.6096834 0.7806332 0.3903166
[51,] 0.5448627 0.9102745 0.4551373
[52,] 0.4527511 0.9055022 0.5472489
[53,] 0.5121330 0.9757340 0.4878670
[54,] 0.4376915 0.8753830 0.5623085
[55,] 0.3807507 0.7615014 0.6192493
[56,] 0.2936381 0.5872763 0.7063619
[57,] 0.5383019 0.9233963 0.4616981
[58,] 0.5701320 0.8597360 0.4298680
[59,] 0.4991608 0.9983215 0.5008392
> postscript(file="/var/www/rcomp/tmp/1c6ds1292182204.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/2ngvd1292182204.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/3qyt01292182204.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/4qyt01292182204.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/5qyt01292182204.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 = 92
Frequency = 1
1 2 3 4 5 6
-2578.99937 -1703.64843 -436.50356 -487.88792 3045.72401 -400.12644
7 8 9 10 11 12
308.05376 2283.76576 -1395.13674 2334.65899 -2054.14288 1631.84263
13 14 15 16 17 18
4692.71309 1569.01241 -3798.42564 -1733.63487 -870.89959 -13.86211
19 20 21 22 23 24
1093.49421 17.53368 912.53853 871.14799 -242.10887 -2174.03206
25 26 27 28 29 30
-2027.19993 -811.40607 -436.24135 2606.50456 -1616.88762 435.11255
31 32 33 34 35 36
76.00776 -3112.33856 1277.63243 -1209.27353 3060.57701 1961.28287
37 38 39 40 41 42
1638.77360 -2775.48000 2200.66212 -2187.52096 -1658.69118 -93.36154
43 44 45 46 47 48
-1049.41110 -1086.90982 239.98177 -2557.87358 2916.84305 -384.97566
49 50 51 52 53 54
-1131.77797 -2149.52956 2008.68462 733.73815 451.55347 35.63504
55 56 57 58 59 60
150.92555 1676.34894 965.03440 -935.94730 1711.50017 2114.67702
61 62 63 64 65 66
-2915.23877 3123.32726 1153.26122 -99.75774 -267.86520 -1094.19996
67 68 69 70 71 72
565.74862 -360.70148 -2650.74311 -313.93725 -2878.67659 -1797.34502
73 74 75 76 77 78
-1714.80362 -304.44692 -1879.86709 -1052.99501 1202.47199 -444.33263
79 80 81 82 83 84
643.83957 -27.48397 650.69273 1811.22468 -2513.99188 -1351.44978
85 86 87 88 89 90
4036.53297 3052.17132 1188.42968 2221.55378 -285.40588 1575.13509
91 92
-1788.65836 609.78544
> postscript(file="/var/www/rcomp/tmp/6jpa31292182204.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 = 92
Frequency = 1
lag(myerror, k = 1) myerror
0 -2578.99937 NA
1 -1703.64843 -2578.99937
2 -436.50356 -1703.64843
3 -487.88792 -436.50356
4 3045.72401 -487.88792
5 -400.12644 3045.72401
6 308.05376 -400.12644
7 2283.76576 308.05376
8 -1395.13674 2283.76576
9 2334.65899 -1395.13674
10 -2054.14288 2334.65899
11 1631.84263 -2054.14288
12 4692.71309 1631.84263
13 1569.01241 4692.71309
14 -3798.42564 1569.01241
15 -1733.63487 -3798.42564
16 -870.89959 -1733.63487
17 -13.86211 -870.89959
18 1093.49421 -13.86211
19 17.53368 1093.49421
20 912.53853 17.53368
21 871.14799 912.53853
22 -242.10887 871.14799
23 -2174.03206 -242.10887
24 -2027.19993 -2174.03206
25 -811.40607 -2027.19993
26 -436.24135 -811.40607
27 2606.50456 -436.24135
28 -1616.88762 2606.50456
29 435.11255 -1616.88762
30 76.00776 435.11255
31 -3112.33856 76.00776
32 1277.63243 -3112.33856
33 -1209.27353 1277.63243
34 3060.57701 -1209.27353
35 1961.28287 3060.57701
36 1638.77360 1961.28287
37 -2775.48000 1638.77360
38 2200.66212 -2775.48000
39 -2187.52096 2200.66212
40 -1658.69118 -2187.52096
41 -93.36154 -1658.69118
42 -1049.41110 -93.36154
43 -1086.90982 -1049.41110
44 239.98177 -1086.90982
45 -2557.87358 239.98177
46 2916.84305 -2557.87358
47 -384.97566 2916.84305
48 -1131.77797 -384.97566
49 -2149.52956 -1131.77797
50 2008.68462 -2149.52956
51 733.73815 2008.68462
52 451.55347 733.73815
53 35.63504 451.55347
54 150.92555 35.63504
55 1676.34894 150.92555
56 965.03440 1676.34894
57 -935.94730 965.03440
58 1711.50017 -935.94730
59 2114.67702 1711.50017
60 -2915.23877 2114.67702
61 3123.32726 -2915.23877
62 1153.26122 3123.32726
63 -99.75774 1153.26122
64 -267.86520 -99.75774
65 -1094.19996 -267.86520
66 565.74862 -1094.19996
67 -360.70148 565.74862
68 -2650.74311 -360.70148
69 -313.93725 -2650.74311
70 -2878.67659 -313.93725
71 -1797.34502 -2878.67659
72 -1714.80362 -1797.34502
73 -304.44692 -1714.80362
74 -1879.86709 -304.44692
75 -1052.99501 -1879.86709
76 1202.47199 -1052.99501
77 -444.33263 1202.47199
78 643.83957 -444.33263
79 -27.48397 643.83957
80 650.69273 -27.48397
81 1811.22468 650.69273
82 -2513.99188 1811.22468
83 -1351.44978 -2513.99188
84 4036.53297 -1351.44978
85 3052.17132 4036.53297
86 1188.42968 3052.17132
87 2221.55378 1188.42968
88 -285.40588 2221.55378
89 1575.13509 -285.40588
90 -1788.65836 1575.13509
91 609.78544 -1788.65836
92 NA 609.78544
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1703.64843 -2578.99937
[2,] -436.50356 -1703.64843
[3,] -487.88792 -436.50356
[4,] 3045.72401 -487.88792
[5,] -400.12644 3045.72401
[6,] 308.05376 -400.12644
[7,] 2283.76576 308.05376
[8,] -1395.13674 2283.76576
[9,] 2334.65899 -1395.13674
[10,] -2054.14288 2334.65899
[11,] 1631.84263 -2054.14288
[12,] 4692.71309 1631.84263
[13,] 1569.01241 4692.71309
[14,] -3798.42564 1569.01241
[15,] -1733.63487 -3798.42564
[16,] -870.89959 -1733.63487
[17,] -13.86211 -870.89959
[18,] 1093.49421 -13.86211
[19,] 17.53368 1093.49421
[20,] 912.53853 17.53368
[21,] 871.14799 912.53853
[22,] -242.10887 871.14799
[23,] -2174.03206 -242.10887
[24,] -2027.19993 -2174.03206
[25,] -811.40607 -2027.19993
[26,] -436.24135 -811.40607
[27,] 2606.50456 -436.24135
[28,] -1616.88762 2606.50456
[29,] 435.11255 -1616.88762
[30,] 76.00776 435.11255
[31,] -3112.33856 76.00776
[32,] 1277.63243 -3112.33856
[33,] -1209.27353 1277.63243
[34,] 3060.57701 -1209.27353
[35,] 1961.28287 3060.57701
[36,] 1638.77360 1961.28287
[37,] -2775.48000 1638.77360
[38,] 2200.66212 -2775.48000
[39,] -2187.52096 2200.66212
[40,] -1658.69118 -2187.52096
[41,] -93.36154 -1658.69118
[42,] -1049.41110 -93.36154
[43,] -1086.90982 -1049.41110
[44,] 239.98177 -1086.90982
[45,] -2557.87358 239.98177
[46,] 2916.84305 -2557.87358
[47,] -384.97566 2916.84305
[48,] -1131.77797 -384.97566
[49,] -2149.52956 -1131.77797
[50,] 2008.68462 -2149.52956
[51,] 733.73815 2008.68462
[52,] 451.55347 733.73815
[53,] 35.63504 451.55347
[54,] 150.92555 35.63504
[55,] 1676.34894 150.92555
[56,] 965.03440 1676.34894
[57,] -935.94730 965.03440
[58,] 1711.50017 -935.94730
[59,] 2114.67702 1711.50017
[60,] -2915.23877 2114.67702
[61,] 3123.32726 -2915.23877
[62,] 1153.26122 3123.32726
[63,] -99.75774 1153.26122
[64,] -267.86520 -99.75774
[65,] -1094.19996 -267.86520
[66,] 565.74862 -1094.19996
[67,] -360.70148 565.74862
[68,] -2650.74311 -360.70148
[69,] -313.93725 -2650.74311
[70,] -2878.67659 -313.93725
[71,] -1797.34502 -2878.67659
[72,] -1714.80362 -1797.34502
[73,] -304.44692 -1714.80362
[74,] -1879.86709 -304.44692
[75,] -1052.99501 -1879.86709
[76,] 1202.47199 -1052.99501
[77,] -444.33263 1202.47199
[78,] 643.83957 -444.33263
[79,] -27.48397 643.83957
[80,] 650.69273 -27.48397
[81,] 1811.22468 650.69273
[82,] -2513.99188 1811.22468
[83,] -1351.44978 -2513.99188
[84,] 4036.53297 -1351.44978
[85,] 3052.17132 4036.53297
[86,] 1188.42968 3052.17132
[87,] 2221.55378 1188.42968
[88,] -285.40588 2221.55378
[89,] 1575.13509 -285.40588
[90,] -1788.65836 1575.13509
[91,] 609.78544 -1788.65836
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1703.64843 -2578.99937
2 -436.50356 -1703.64843
3 -487.88792 -436.50356
4 3045.72401 -487.88792
5 -400.12644 3045.72401
6 308.05376 -400.12644
7 2283.76576 308.05376
8 -1395.13674 2283.76576
9 2334.65899 -1395.13674
10 -2054.14288 2334.65899
11 1631.84263 -2054.14288
12 4692.71309 1631.84263
13 1569.01241 4692.71309
14 -3798.42564 1569.01241
15 -1733.63487 -3798.42564
16 -870.89959 -1733.63487
17 -13.86211 -870.89959
18 1093.49421 -13.86211
19 17.53368 1093.49421
20 912.53853 17.53368
21 871.14799 912.53853
22 -242.10887 871.14799
23 -2174.03206 -242.10887
24 -2027.19993 -2174.03206
25 -811.40607 -2027.19993
26 -436.24135 -811.40607
27 2606.50456 -436.24135
28 -1616.88762 2606.50456
29 435.11255 -1616.88762
30 76.00776 435.11255
31 -3112.33856 76.00776
32 1277.63243 -3112.33856
33 -1209.27353 1277.63243
34 3060.57701 -1209.27353
35 1961.28287 3060.57701
36 1638.77360 1961.28287
37 -2775.48000 1638.77360
38 2200.66212 -2775.48000
39 -2187.52096 2200.66212
40 -1658.69118 -2187.52096
41 -93.36154 -1658.69118
42 -1049.41110 -93.36154
43 -1086.90982 -1049.41110
44 239.98177 -1086.90982
45 -2557.87358 239.98177
46 2916.84305 -2557.87358
47 -384.97566 2916.84305
48 -1131.77797 -384.97566
49 -2149.52956 -1131.77797
50 2008.68462 -2149.52956
51 733.73815 2008.68462
52 451.55347 733.73815
53 35.63504 451.55347
54 150.92555 35.63504
55 1676.34894 150.92555
56 965.03440 1676.34894
57 -935.94730 965.03440
58 1711.50017 -935.94730
59 2114.67702 1711.50017
60 -2915.23877 2114.67702
61 3123.32726 -2915.23877
62 1153.26122 3123.32726
63 -99.75774 1153.26122
64 -267.86520 -99.75774
65 -1094.19996 -267.86520
66 565.74862 -1094.19996
67 -360.70148 565.74862
68 -2650.74311 -360.70148
69 -313.93725 -2650.74311
70 -2878.67659 -313.93725
71 -1797.34502 -2878.67659
72 -1714.80362 -1797.34502
73 -304.44692 -1714.80362
74 -1879.86709 -304.44692
75 -1052.99501 -1879.86709
76 1202.47199 -1052.99501
77 -444.33263 1202.47199
78 643.83957 -444.33263
79 -27.48397 643.83957
80 650.69273 -27.48397
81 1811.22468 650.69273
82 -2513.99188 1811.22468
83 -1351.44978 -2513.99188
84 4036.53297 -1351.44978
85 3052.17132 4036.53297
86 1188.42968 3052.17132
87 2221.55378 1188.42968
88 -285.40588 2221.55378
89 1575.13509 -285.40588
90 -1788.65836 1575.13509
91 609.78544 -1788.65836
> 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/rcomp/tmp/7ch961292182204.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/84qr91292182204.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/94qr91292182204.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/rcomp/tmp/104qr91292182204.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/www/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/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/rcomp/tmp/11fsmr1292182204.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/rcomp/tmp/12tj2i1292182204.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/rcomp/tmp/13augw1292182204.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/rcomp/tmp/14wcxk1292182204.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/rcomp/tmp/15zddq1292182204.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/rcomp/tmp/162duw1292182204.tab")
+ }
>
> try(system("convert tmp/1c6ds1292182204.ps tmp/1c6ds1292182204.png",intern=TRUE))
character(0)
> try(system("convert tmp/2ngvd1292182204.ps tmp/2ngvd1292182204.png",intern=TRUE))
character(0)
> try(system("convert tmp/3qyt01292182204.ps tmp/3qyt01292182204.png",intern=TRUE))
character(0)
> try(system("convert tmp/4qyt01292182204.ps tmp/4qyt01292182204.png",intern=TRUE))
character(0)
> try(system("convert tmp/5qyt01292182204.ps tmp/5qyt01292182204.png",intern=TRUE))
character(0)
> try(system("convert tmp/6jpa31292182204.ps tmp/6jpa31292182204.png",intern=TRUE))
character(0)
> try(system("convert tmp/7ch961292182204.ps tmp/7ch961292182204.png",intern=TRUE))
character(0)
> try(system("convert tmp/84qr91292182204.ps tmp/84qr91292182204.png",intern=TRUE))
character(0)
> try(system("convert tmp/94qr91292182204.ps tmp/94qr91292182204.png",intern=TRUE))
character(0)
> try(system("convert tmp/104qr91292182204.ps tmp/104qr91292182204.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.470 1.660 5.122