wsMean.R
In FitARMA: Fit ARMA or ARIMA Using Fast MLE Algorithm

"cleanws" <-
function(){
rm(list=ls(envir=.GlobalEnv), envir=.GlobalEnv)
load(file=".Rdata", envir=.GlobalEnv)
}

"er" <-
structure(c(1.01356909276102, 1.03398545497203, 1.03641136673971, 
1.02649357536719, 1.02001438440014, 0.909797131547903, 0.698834269339134, 
0.992303506072454, 1.03289855443273, 1.04721300604513, 1.04093797971543, 
1.03279016815024, 0.915128727361327, 0.704032758694548, 0.993542004521725, 
0.976320256546273, 1.01459655415719, 1.04093748845283, 1.07439390725605, 
0.965380426882639, 0.740927754685409, 0.771138511387702, 1.05025835186947, 
1.01341805204479, 1.02668114916583, 1.04053546121134, 1.08965097336703, 
0.913403706430706, 0.311182734393226, 0.693437255256049, 0.881806921945306, 
1.03212636815038, 0.986675278567097, 0.966766110229769, 0.782671724640663, 
0.00822672942647032, 0.0523315847827543, 0.181472474458157, 1.03212637683340, 
0.458251574415404, 1.14521370492097, 1.02957619545235, 0.00138762682955874, 
0.0100676795452076, 0.0485714375517598, 1.03212982018543, 0.216725361714229, 
0.905718582472695, 1.11207998204141), .Dim = as.integer(c(7, 
7)), .Dimnames = list(c("phi=-0.975", "phi=-0.9", "phi=-0.5", 
"phi=0", "phi=0.5", "phi=0.9", "phi=0.975"), c("theta=-0.975", 
"theta=-0.9", "theta=-0.5", "theta=0", "theta=0.5", "theta=0.9", 
"theta=0.975")))
"erarima" <-
structure(c(1.01003730213983, 0.973412787845336, 1.03765273766072, 
0.99105013396048, 1.02377789319720, 0.908686324499887, 0.697092592331931, 
1.04599850246196, 1.03945680295961, 1.04807382519451, 0.99105013396048, 
1.03701487763842, 0.917099134575337, 0.705034184587719, 0.958623836004144, 
0.97633463927683, 1.02438827739318, 0.99105013396048, 1.07440894409782, 
0.965391759352532, 0.741002750567826, 0.850943796459182, 1.05052572514667, 
1.01346200457931, 1.02816716714078, 1.04054023962569, 1.08975106818921, 
0.913336098953355, 0.31127074954211, 0.693434564200314, 0.881816345743621, 
0.99105013396048, 1.05235238300339, 0.966649990619832, 0.783057480066306, 
0.00592165529941912, 0.0376607531651018, 0.136708771856345, 0.99105013396048, 
0.400422946892324, 1.04986227193339, 1.03656099831922, 0.000639489335717399, 
0.00458406843763869, 0.0213439310678546, 0.99105013396048, 0.173372676019292, 
0.929407214866703, 1.10900661221307), .Dim = as.integer(c(7, 
7)), .Dimnames = list(c("phi=-0.975", "phi=-0.9", "phi=-0.5", 
"phi=0", "phi=0.5", "phi=0.9", "phi=0.975"), c("theta=-0.975", 
"theta=-0.9", "theta=-0.5", "theta=0", "theta=0.5", "theta=0.9", 
"theta=0.975")))
"GetARMAMeanFit" <-
function(p,q,z,pApprox=30,init=0){
mz<-mean(z)
y<-z-mz
mu<-0
if (q>0)
   pApp <- pApprox
else
   pApp <- p
ans<-GetFitARMA(y,p,q,init=init)
if (ans$convergence != 0) warning("GetARMAMeanFit: convergence = ", ans$convergence)
pinit<-c(ARToPacf(ans$phi),ARToPacf(ans$theta))
if (max(abs(pinit))>0.99){
    message("GetARMAMeanFit: bailing out. Used sample mean, pinit=", pinit)
    ZZ<<-y #save problem data
    stop()
    return(c(ans, mean=mz, NumIterations=0, etol=NULL))
    }
logL<-ans$loglikelihood
g<-TacvfARMA(ans$phi,ans$theta,pApp)
coefAR<-PacfDL(g, LinearPredictor=T)$ARCoefficients
mu<-GetARMeanMLE(y, coefAR)
x<-y-mu
ans<-GetFitARMA(x,p,q,init=pinit)
etol<-abs(ans$loglikelihood-logL)
iter<-1
while(etol>1e-3 && iter<5){
    iter<-iter+1
    logL<-ans$loglikelihood
    g<-TacvfARMA(ans$phi,ans$theta,pApp)
    coefAR<-PacfDL(g, LinearPredictor=T)$ARCoefficients
    mu<-GetARMeanMLE(y, coefAR)
    x<-y-mu
    pinit<-c(ARToPacf(ans$phi),ARToPacf(ans$theta))
    ans<-GetFitARMA(x,p,q,init=pinit)    
    etol<-abs(ans$loglikelihood-logL)
    }
c(ans, mean=mz+mu, NumIterations=iter, etol=etol)
}

"last.warning" <-
structure(list("NAs introduced by coercion" = NULL), .Names = "NAs introduced by coercion")
"loadws" <-
function(name, d="2007"){
name<-paste(d,name,sep="/")
ws <- paste(paste("d:/r",name,sep="/"),".Rdata",sep="/")
wsH <- paste(paste("d:/r",name,sep="/"),".RHistory",sep="/")
cat(ws, fill=T)
load(ws, .GlobalEnv)
loadhistory(wsH)
cat(paste("loaded: ", ws), fill=T)
cat(paste("loaded: ", wsH), fill=T)
}

"MeanEffSim" <-
function(phi,theta,n,NSIM){
mout<-matrix(numeric(NSIM*3), ncol=3, nrow=NSIM)
i<-0
set.seed(start.seed)
while (i < NSIM) {
    z <- SimulateGaussianARMA(phi,theta,n)
    zm <- mean(z)
    #
    #arima
    zmarima <- tryCatch(coef(arima(z, c(1,0,1)))[3],error = function(e) NA, NA)
    #
    #FitARMA
    zmmle <- GetARMAMeanFit(1,1,z,init=c(ARToPacf(phi),ARToPacf(theta)))$mean
    out<-c(zm,zmmle,zmarima)
    if (length(out)!=3) message("zmmle=",zmmle, ", zmarima=", zmarima)
    i<-i+1
    mout[i,] <- out
   }
   mout
}

"RunMeanEffSim" <-
function(NSIM){
er<<-0
erarima<<-0
pars<-c(-0.975, -0.9,-0.5,0,0.5,0.9, 0.975)
np<-length(pars)
er<-erarima<-matrix(numeric(np*np),nrow=np,ncol=np)
for (i in 1:np){
    for (j in 1:np){
#there are some convergence problems with the built-in R function arima
        out<-MeanEffSim(pars[i],pars[j],100,NSIM)
        v<-apply(out^2,MARGIN=2,FUN=function(x) mean(x, na.rm=T))
        er[i,j]<- v[2]/v[1]
        erarima[i,j]<- v[3]/v[1]
        }
    }
rn<-paste("phi=",pars,sep="")
cn<-paste("theta=",pars,sep="")
dimnames(er)<-list(rn,cn)
dimnames(erarima)<-list(rn,cn)
er<<-er
erarima<<-erarima
NULL
}

"savews" <-
function(name=.WSID, d="2007"){
if (!exists(".WSID") && !is.character(name)) stop(".WSID not defined!")
.WSID<<-name
name<-paste(d,name,sep="/")
wsRdata <<- paste(paste("d:/r", name, sep="/"), ".Rdata", sep="/")
wsR <<- paste(paste("d:/r", name, sep="/"), "ws.R", sep="/")
save.image(wsRdata)
dump(ls(envir=.GlobalEnv),wsR)
wsRHistory <<- paste(paste("d:/r", name, sep="/"), ".RHistory", sep="/")
savehistory(wsRHistory)
cat(paste("saved: ", wsRdata), fill=T)
cat(paste("saved: ", wsR), fill=T)
cat(paste("saved: ", wsRHistory), fill=T)
}

"start.seed" <-
as.integer(c(403, 624, 1225564623, -987490876, -926653963, 540074546, 
617851915, -381135920, -822301295, 1119377950, 609333255, -345943908, 
-1096714003, 1398726346, 1930602947, -1034692568, -1220141047, 
1926306870, -1407691457, 1628403316, -119947547, 334790242, -476696965, 
236615552, 448999809, -1991036594, 1651339639, -773986484, 985469917, 
-1220009734, -1888937421, 1202819544, -9321991, 384498022, 1111090351, 
-576191708, 181884885, -172218734, 2083671275, 1327138608, 944484721, 
-1643064194, 1120046823, -436917764, -1032820019, -734073046, 
397715107, -826099832, 751230953, -729210218, 1260933151, -1955021868, 
-1928581947, -1110779198, 592381787, 1359184608, -1873466015, 
350503854, -1779955625, -611182420, 1474985405, -857323174, 259804947, 
94521656, 165968345, -92218, -2074437745, 768385156, -1316536907, 
1159961330, -950837813, 928017040, -882989743, 1317043934, -499851833, 
-1319128228, -1726329683, 1017196426, -335080573, 1898661608, 
422154185, -765568778, -1727546625, -1431391948, 1041730213, 
1971939106, -2140778437, 1413234240, -1046013631, -1470593522, 
-742387913, 1623783948, -1347496035, 1737662906, 73135091, 486060184, 
-2068504135, -1719966170, -1842955665, -1159074332, 2085360533, 
2038382418, 7265963, -658375184, 1798146353, -1398843074, -1402954585, 
-1913066308, 1092034189, 1788714986, 21098595, 1264944712, 241577897, 
-404177066, 1181652447, -1710664044, 695457925, -330816638, 51644699, 
-2070107744, -810486495, 1112012910, -1280473577, 937069420, 
1585586557, 1801792026, 1296042195, 629983224, -2377319, -990288762, 
-900313777, -1235752124, 1721620853, 166138802, -1111699573, 
1447510352, -156213999, -586849378, -1493316729, 1746458140, 
1960763501, -934526902, -2070069949, 1799976360, 443858825, -591374922, 
-535502657, 1625336820, -1466362267, -551612446, 1334850043, 
989644032, -680869631, -1387619634, 694709495, -484520756, 1021074269, 
1215685242, -446657101, 590778200, -2024783495, -1341091098, 
1738625071, 2104400036, -1477144747, 2037585938, 829763691, -1095224144, 
1222582513, -800416258, 868510311, -669552772, -1427533235, 1271206058, 
-1085212125, 1252988168, -951238807, -998433770, -827155553, 
848081236, 1358906437, 388377666, -1508717861, -1137405856, -238256927, 
2146987310, -2069314601, -1863467476, -550140611, -86193446, 
-470449517, -2050095432, -1559578279, -732368570, -2134916337, 
-1719273980, -1278724811, 1463503986, 701497675, 216848400, 951178449, 
1124534366, 275544391, 575453404, 363803693, 1998590218, 378873603, 
-814848920, 451391305, -21556618, 1459603071, 1852139188, -199335387, 
1769306274, -489434181, 958139328, 869148865, 485060494, 1898351287, 
463429516, -1778016483, 179430202, 2083972979, 690695704, 1459822905, 
-144015450, 141141487, -1074396316, 965790485, 1045934290, 285557291, 
666708848, -1725925199, -1110269250, 1454241831, 923841084, -1549285875, 
1484079466, 195150819, 1264542664, -1757671639, 923154646, -1716231841, 
-1614623724, -1620152315, -1222314750, 2114669723, -829737184, 
-1368931167, -1296719378, -157695593, 160149740, 1841604861, 
-1945293926, -44254125, 1428555128, 590444825, 719142406, -849938225, 
-770260796, 616976629, -706722510, -313724149, -503238960, 958598289, 
-1790612194, 2084628231, -1296344164, -11843603, -1982029366, 
1001312451, 2086277928, 977428233, 1694666550, -1819775937, 2113724788, 
-1670943259, -313745054, -1947126405, -2057694592, 1855017089, 
1005913166, 2025482359, -1928678836, 787129053, 520528890, -761331405, 
-1114207016, 59621625, -867619738, 2050997167, -582996444, -1662990635, 
-789771886, 1625232363, -91168208, -474902415, -369662082, 1390233063, 
-576408324, -1894663731, 829305386, 1609846179, -1942588792, 
2113052393, -867953258, 489940767, -328489260, 1887525829, 84180418, 
-1128427941, 1382063584, 2101529697, 2134873774, -1220509865, 
-2072747092, 1575980221, -273265574, -2093631981, -2108373960, 
1880094937, -2058990906, -1780748657, 449465220, 89664693, -292160014, 
-1443650357, 474236304, 1606681681, -1467969058, 138089671, -1406878116, 
1970477997, 21099146, 1303001731, 281838056, 1480904393, -60634122, 
-345058817, -93905868, -583779931, 56920610, 1554536251, 376887616, 
-546031551, 254630158, -840694217, 2048967948, 994797213, -1238096710, 
-1302799629, 792243096, 1555046585, 1265408294, -2099015313, 
-267577116, -39550315, -101506478, -1564301909, -571254544, 1884448817, 
-1988563906, 638735271, -1097631812, -1863881331, 930079466, 
-167209117, 200556872, 993062569, -789138858, -1896796961, -481869420, 
-537393275, -108738942, 1407815707, -1636514656, -1776446431, 
1247169390, 878509335, -1536698772, -1015664515, -1231540966, 
524316627, -1139129608, 963000473, 1616123782, -2146524081, -295659964, 
1631438965, -904104270, -1048308085, -1032446896, -732644335, 
293107358, -1843723641, 1632166172, 2078715757, -1643117750, 
1715954755, -343560024, 347012745, 1905772726, 1953691583, 241440500, 
-1304115867, 75313890, 649672955, -1556979712, -1640570879, 1532128718, 
-1100782601, -442394676, -1360532899, -1057331846, -1524337485, 
-1932424616, -232112007, 1368704486, -1475008721, -753089628, 
1201405525, 1110047506, 409991019, 950308784, 1187184625, -1860751106, 
-1911742105, -1935869316, -1930895027, -1959111766, -945905373, 
-1990668280, 1320615529, 1373507350, -458476897, 253074516, -913149115, 
1298517826, -179351077, -913855648, -256369695, 951697454, -1683014953, 
-869922516, 1914213437, 829607386, 993847699, 1899397560, -285984679, 
-121199546, -240182769, -2019974908, 70722613, 1362341746, 1438533707, 
-1688816880, 2123707345, 619518814, -1214205881, -454573092, 
-697957587, -519646198, 1598443011, 725983080, -817828279, 828474742, 
72670591, -1531719244, -682028763, 156670882, 2078530235, -1572034880, 
-2103883839, -1524350322, 1475903927, -1840436596, 1031810589, 
-266770886, -188625293, -1524553448, 331096121, 2072097446, 798260463, 
556740196, 646396437, -229534766, -1012463317, 828671600, 784553905, 
-1248709186, 143503143, -1165935300, 651564301, 227378282, -1904842013, 
-2094784824, -129708503, 465185750, -735774625, -1164527852, 
-921657595, 2036832258, 313447323, -1436986848, 1154638753, 791278830, 
-621332329, 510589932, -40454147, 1896230554, 15410003, -798392200, 
-1065748455, 1264447750, 276647887, -516592700, 2047098869, 748398642, 
1214396939, 765856208, 105212817, -140607458, -695460345, 163669660, 
145823469, 182871242, -765003837, -1362342360, -1480942071, 1753219638, 
819233599, 1746291828, -936305435, -357514142, -1377289093, 1250274688, 
609972097, 828561230, 1751342967, 48463180, 1523855837, -1269750022, 
-1427052493, 115836888, -786055175, 636706662, 522293935, 889477412, 
103167445, 321514642, 1713887979, -1359790800, -1340903567, 1806301822, 
-749470489, 2063613948, -1032910643, 1596552490, -1001392989, 
940977544, 932863465, -1152710506, -753172961, -259354156, 976390853, 
-1236655934, -574089893, -676742944, 170683233, -765007442, -1611336105, 
2114105004, -1179608131, 1175605082, 1510677779, -931971272, 
-1648920615, 349830598, -1036434033, -1242302844, 241507253, 
-988520206, 992996299, -1071374192, -752524463, 1582081246, 211635143, 
1653983580, 1451748013, 477017482, 426336643, 349813992, 2111573449, 
62078710, 1337059583, -1018460364, -1064507227, 1095478562, -830054853, 
-1835174848, -616736959, 80620558, 2103704887, 2049216524, 983823773, 
1046587322, -1854815757, 14984856, -98099271, 1839844390, 1014786159, 
693912548, 305721749, 1836254546, -2119013205, 1477489648, 309544753, 
-400654530, -333444441, -1059454276, -2089567091, -523360790, 
-1561641373, -1494287288, -864571991, 2102437206))
"wsR" <-
"d:/r/2007/FitARMA/Mean/ws.R"
"wsRdata" <-
"d:/r/2007/FitARMA/Mean/.Rdata"
"wsRHistory" <-
"d:/r/2007/FitARMA/Mean/.RHistory"