R : Copyright 2005, The R Foundation for Statistical Computing Version 2.1.1 (2005-06-20), 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 a HTML browser interface to help. Type 'q()' to quit R. > ### *

> ### > attach(NULL, name = "CheckExEnv") > assign(".CheckExEnv", as.environment(2), pos = length(search())) # base > ## add some hooks to label plot pages for base and grid graphics > setHook("plot.new", ".newplot.hook") > setHook("persp", ".newplot.hook") > setHook("grid.newpage", ".gridplot.hook") > > assign("cleanEx", + function(env = .GlobalEnv) { + rm(list = ls(envir = env, all.names = TRUE), envir = env) + RNGkind("default", "default") + set.seed(1) + options(warn = 1) + delayedAssign("T", stop("T used instead of TRUE"), + assign.env = .CheckExEnv) + delayedAssign("F", stop("F used instead of FALSE"), + assign.env = .CheckExEnv) + sch <- search() + newitems <- sch[! sch %in% .oldSearch] + for(item in rev(newitems)) + eval(substitute(detach(item), list(item=item))) + missitems <- .oldSearch[! .oldSearch %in% sch] + if(length(missitems)) + warning("items ", paste(missitems, collapse=", "), + " have been removed from the search path") + }, + env = .CheckExEnv) > assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now > assign("ptime", proc.time(), env = .CheckExEnv) > grDevices::postscript("micEcon-Examples.ps") > assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv) > options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly")) > options(warn = 1) > library('micEcon') Loading required package: systemfit This package is still an alpha version. This means that this package has not been thoroughly tested yet. Thus, there might be several bugs in it and I do not guarantee that all results are correct. Please send bug reports to ahenningsen@agric-econ.uni-kiel.de > > assign(".oldSearch", search(), env = .CheckExEnv) > assign(".oldNS", loadedNamespaces(), env = .CheckExEnv) > cleanEx(); ..nameEx <- "aidsBestA0" > > ### * aidsBestA0 > > flush(stderr()); flush(stdout()) > > ### Name: aidsBestA0 > ### Title: Find 'best' Value for Coefficient alpha 0 > ### Aliases: aidsBestA0 > ### Keywords: models > > ### ** Examples > > data( Blanciforti86 ) > bestA0 <- aidsBestA0( c( "pFood1", "pFood2", "pFood3", "pFood4" ), + c( "wFood1", "wFood2", "wFood3", "wFood4" ), "xFood", + data = Blanciforti86, method = "MK:L" ) # may take some time > bestA0$alpha0 [1] -1.041667 > plot( bestA0$allValues ) # this should be convex > > > > cleanEx(); ..nameEx <- "aidsCalc" > > ### * aidsCalc > > flush(stderr()); flush(stdout()) > > ### Name: aidsCalc > ### Title: Shares and Quantities of the Almost Ideal Demand System > ### Aliases: aidsCalc > ### Keywords: models > > ### ** Examples > > data( Blanciforti86 ) > pNames <- c( "pFood1", "pFood2", "pFood3", "pFood4" ) > wNames <- c( "wFood1", "wFood2", "wFood3", "wFood4" ) > > ## LA-AIDS > estResult <- aidsEst( pNames, wNames, "xFood", + data = Blanciforti86, method = "LA:L" ) > > lnp <- aidsPx( "L", pNames, wNames, Blanciforti86 ) > > fitted <- aidsCalc( pNames, "xFood", Blanciforti86, + coef = estResult$coef, lnp = lnp ) > > fitted$shares # equal to estResult$wFitted w1 w2 w3 w4 1947 0.3164903 0.1879093 0.1327508 0.3628496 1948 0.3162674 0.1716257 0.1351846 0.3769224 1949 0.3069036 0.1828754 0.1358315 0.3743895 1950 0.3124188 0.1687368 0.1367450 0.3820994 1951 0.3141931 0.1641911 0.1369036 0.3847122 1952 0.3026762 0.1813374 0.1368023 0.3791841 1953 0.3059778 0.1841025 0.1367688 0.3731510 1954 0.2991363 0.1828327 0.1390867 0.3789442 1955 0.2996730 0.1957473 0.1377669 0.3668128 1956 0.2937063 0.2106971 0.1373492 0.3582474 1957 0.2982332 0.1977530 0.1379634 0.3660504 1958 0.2934549 0.1955206 0.1382332 0.3727912 1959 0.3007079 0.2002430 0.1364914 0.3625576 1960 0.3001759 0.2063471 0.1360103 0.3574667 1961 0.2986112 0.2058844 0.1367369 0.3587675 1962 0.3038012 0.2044625 0.1357573 0.3559789 1963 0.2978529 0.2159930 0.1354618 0.3506924 1964 0.2943856 0.2239734 0.1351494 0.3464916 1965 0.3033082 0.2151138 0.1340139 0.3475641 1966 0.3128774 0.2061202 0.1328848 0.3481176 1967 0.3179425 0.2099640 0.1318884 0.3402051 1968 0.3146346 0.2190852 0.1307915 0.3354888 1969 0.3235171 0.2082345 0.1300144 0.3382340 1970 0.3274147 0.2059551 0.1296598 0.3369704 1971 0.3188726 0.2122821 0.1310292 0.3378161 1972 0.3258681 0.2076161 0.1290387 0.3374771 1973 0.3321180 0.1951422 0.1279959 0.3447438 1974 0.3178235 0.2134358 0.1313599 0.3373808 1975 0.3231290 0.2041811 0.1321924 0.3404975 1976 0.3295624 0.2055246 0.1303644 0.3345486 1977 0.3116443 0.2164798 0.1329089 0.3389669 1978 0.3175834 0.2116033 0.1313051 0.3395082 > fitted$quant # equal to estResult$qFitted q1 q2 q3 q4 1947 1.649878 1.0837220 0.7925990 1.941456 1948 1.525730 1.0095206 0.7680357 1.980181 1949 1.510062 1.0257768 0.7511418 2.015911 1950 1.469941 0.9852042 0.7417491 2.014511 1951 1.462260 0.9714752 0.7507000 1.984172 1952 1.459167 1.0022643 0.7531800 2.014083 1953 1.559751 1.0503940 0.7381119 1.910492 1954 1.514191 1.0340711 0.7147292 1.868851 1955 1.629848 1.1021359 0.7060117 1.842886 1956 1.700461 1.1493028 0.7048193 1.801643 1957 1.628892 1.1061731 0.6964051 1.857829 1958 1.505459 1.0509217 0.7043415 1.986089 1959 1.620907 1.1053370 0.6968070 1.975804 1960 1.674387 1.1344816 0.6963085 1.946532 1961 1.673472 1.1350675 0.6898018 1.912904 1962 1.698802 1.1441329 0.6846488 1.962909 1963 1.729512 1.1678974 0.6891085 1.944480 1964 1.769383 1.1890470 0.6984280 1.899020 1965 1.762608 1.1774141 0.7022441 1.969255 1966 1.777030 1.1721185 0.7104917 1.998848 1967 1.902850 1.2267437 0.7067480 1.928598 1968 1.927684 1.2403213 0.7306295 1.924628 1969 1.920550 1.2251279 0.7422468 1.943065 1970 1.962722 1.2385345 0.7451289 1.907831 1971 1.943904 1.2399702 0.7340388 1.872732 1972 1.912846 1.2187063 0.7574571 1.980991 1973 1.777965 1.1482955 0.7715024 2.179717 1974 1.900267 1.2314155 0.6944018 2.007542 1975 1.918023 1.2319442 0.6776553 1.961099 1976 2.044207 1.2742196 0.7198098 1.833512 1977 2.002046 1.2651477 0.7437538 1.673699 1978 1.941931 1.2357669 0.7496298 1.817974 > > ## AIDS > estResult <- aidsEst( pNames, wNames, "xFood", + data = Blanciforti86, method = "MK:L" ) > > fitted <- aidsCalc( pNames, "xFood", Blanciforti86, + coef = estResult$coef ) > > fitted$shares # equal to estResult$wFitted w1 w2 w3 w4 1947 0.3168496 0.1881344 0.1327943 0.3622217 1948 0.3166410 0.1718889 0.1352359 0.3762342 1949 0.3072543 0.1831310 0.1358084 0.3738063 1950 0.3127573 0.1689913 0.1367687 0.3814827 1951 0.3143707 0.1644348 0.1369899 0.3842045 1952 0.3030071 0.1816287 0.1367706 0.3785936 1953 0.3060261 0.1842270 0.1367863 0.3729606 1954 0.2990745 0.1829126 0.1390726 0.3789402 1955 0.2997220 0.1957476 0.1376886 0.3668418 1956 0.2936027 0.2106254 0.1372437 0.3585282 1957 0.2981594 0.1977223 0.1378941 0.3662241 1958 0.2930681 0.1955809 0.1382440 0.3731070 1959 0.3003173 0.2002106 0.1365153 0.3629568 1960 0.2996244 0.2062529 0.1360538 0.3580689 1961 0.2981354 0.2057773 0.1367455 0.3593418 1962 0.3033585 0.2043474 0.1357827 0.3565114 1963 0.2968508 0.2157926 0.1355732 0.3517834 1964 0.2930601 0.2237156 0.1353122 0.3479122 1965 0.3024397 0.2149473 0.1341375 0.3484755 1966 0.3126016 0.2060478 0.1329413 0.3484094 1967 0.3178216 0.2098169 0.1319116 0.3404499 1968 0.3140480 0.2189144 0.1309187 0.3361189 1969 0.3235480 0.2081741 0.1300718 0.3382061 1970 0.3277163 0.2059033 0.1296747 0.3367057 1971 0.3189396 0.2121799 0.1310326 0.3378479 1972 0.3259463 0.2076069 0.1291171 0.3373296 1973 0.3335990 0.1954761 0.1278282 0.3430966 1974 0.3173425 0.2132394 0.1314544 0.3379637 1975 0.3231957 0.2039950 0.1321749 0.3406344 1976 0.3304438 0.2054372 0.1302159 0.3339030 1977 0.3134837 0.2165522 0.1324391 0.3375250 1978 0.3179561 0.2115582 0.1312430 0.3392428 > fitted$quant # equal to estResult$qFitted q1 q2 q3 q4 1947 1.651751 1.0850202 0.7928585 1.938096 1948 1.527533 1.0110687 0.7683273 1.976566 1949 1.511787 1.0272105 0.7510140 2.012771 1950 1.471534 0.9866900 0.7418775 2.011259 1951 1.463087 0.9729176 0.7511733 1.981554 1952 1.460763 1.0038740 0.7530054 2.010947 1953 1.559998 1.0511045 0.7382064 1.909518 1954 1.513878 1.0345229 0.7146568 1.868831 1955 1.630115 1.1021373 0.7056106 1.843031 1956 1.699861 1.1489118 0.7042776 1.803055 1957 1.628489 1.1060017 0.6960553 1.858711 1958 1.503475 1.0512459 0.7043965 1.987771 1959 1.618801 1.1051579 0.6969289 1.977980 1960 1.671311 1.1339637 0.6965310 1.949811 1961 1.670806 1.1344767 0.6898449 1.915967 1962 1.696327 1.1434885 0.6847770 1.965844 1963 1.723694 1.1668140 0.6896753 1.950529 1964 1.761416 1.1876782 0.6992690 1.906806 1965 1.757561 1.1765030 0.7028921 1.974419 1966 1.775463 1.1717064 0.7107938 2.000523 1967 1.902127 1.2258841 0.7068723 1.929986 1968 1.924091 1.2393546 0.7313404 1.928243 1969 1.920733 1.2247727 0.7425745 1.942905 1970 1.964531 1.2382231 0.7452145 1.906332 1971 1.944312 1.2393732 0.7340578 1.872908 1972 1.913305 1.2186528 0.7579175 1.980125 1973 1.785894 1.1502602 0.7704915 2.169302 1974 1.897391 1.2302824 0.6949012 2.011011 1975 1.918418 1.2308213 0.6775656 1.961888 1976 2.049674 1.2736779 0.7189901 1.829973 1977 2.013863 1.2655707 0.7411247 1.666579 1978 1.944210 1.2355038 0.7492749 1.816552 > > > > > cleanEx(); ..nameEx <- "aidsEla" > > ### * aidsEla > > flush(stderr()); flush(stdout()) > > ### Name: aidsEla > ### Title: Elasticities of the AIDS model > ### Aliases: aidsEla > ### Keywords: models > > ### ** Examples > > data( Blanciforti86 ) > estResult <- aidsEst( c( "pFood1", "pFood2", "pFood3", "pFood4" ), + c( "wFood1", "wFood2", "wFood3", "wFood4" ), "xFood", + data = Blanciforti86, method = "LA:L" ) > wMeans <- colMeans( Blanciforti86[ , c( "wFood1", "wFood2", + "wFood3", "wFood4" ) ] ) > aidsEla( estResult$coef, wMeans, formula = "Ch" ) $formula [1] "Ch" $exp [1] 2.0521948 1.2339508 0.4085361 0.1720357 $hicks pFood1 pFood2 pFood3 pFood4 wFood1 -0.3935093 -0.273672746 0.1048837 0.5622983 wFood2 -0.4239348 0.002774164 0.1525367 0.2686240 wFood3 0.2426583 0.227820959 -0.7599730 0.2894937 wFood4 0.4913202 0.151521804 0.1093327 -0.7521747 $marshall pFood1 pFood2 pFood3 pFood4 wFood1 -1.0303927 -0.6848152 -0.17039515 -0.1665918 wFood2 -0.8068822 -0.2444390 -0.01298398 -0.1696455 wFood3 0.1158721 0.1459737 -0.81477349 0.1443915 wFood4 0.4379302 0.1170557 0.08625602 -0.8132776 > > > > cleanEx(); ..nameEx <- "aidsEst" > > ### * aidsEst > > flush(stderr()); flush(stdout()) > > ### Name: aidsEst > ### Title: Estimation of the Almost Ideal Demand System (AIDS) > ### Aliases: aidsEst > ### Keywords: models > > ### ** Examples > > ## Repeating the demand analysis of Blanciforti, Green & King (1986) > data( Blanciforti86 ) > estResult <- aidsEst( c( "pFood1", "pFood2", "pFood3", "pFood4" ), + c( "wFood1", "wFood2", "wFood3", "wFood4" ), "xFood", + data = Blanciforti86, method = "LA:SL", elaFormula = "Ch", + maxiter = 1, rcovformula = 1, tol = 1e-7 ) > print( estResult ) Demand analysis with the Almost Ideal Demand System (AIDS) Estimation Method: Linear Approximation (LA) with lagged Stone Index (SL) Estimated Coefficients alpha wFood1 wFood2 wFood3 wFood4 -0.2529649 0.1177580 0.2596590 0.8755479 beta wFood1 wFood2 wFood3 wFood4 0.32757754 0.05134219 -0.07384703 -0.30507269 gamma pFood1 pFood2 pFood3 pFood4 wFood1 0.10775725 -0.140800303 -0.0106393511 0.0436824073 wFood2 -0.14080030 0.160963783 -0.0054313908 -0.0147320885 wFood3 -0.01063935 -0.005431391 0.0162411167 -0.0001703748 wFood4 0.04368241 -0.014732088 -0.0001703748 -0.0287799440 Demand Elasticities (Formula of Chalfant / Goddard) Expenditure Elasticities wFood1 wFood2 wFood3 wFood4 2.0541727 1.2550935 0.4495368 0.1378080 Marshallian (uncompensated) Price Elasticities pFood1 pFood2 pFood3 pFood4 wFood1 -0.98080543 -0.66527886 -0.17566022 -0.2324282 wFood2 -0.77883457 -0.25159413 -0.06120775 -0.1634570 wFood3 0.09174618 0.07030455 -0.80509006 0.1935025 wFood4 0.39137531 0.13189614 0.11518530 -0.7762648 Hicksian (compensated) Price Elasticities pFood1 pFood2 pFood3 pFood4 wFood1 -0.3424842 -0.251839394 0.09991601 0.4944076 wFood2 -0.3888222 0.001016173 0.10716851 0.2806375 wFood3 0.2314369 0.160781984 -0.74478273 0.3525638 wFood4 0.4341983 0.159632499 0.13367285 -0.7275036 > > ## Repeating the evaluation of different elasticity formulas of > ## Green & Alston (1990) > data( Blanciforti86 ) > pNames <- c( "pFood1", "pFood2", "pFood3", "pFood4" ) > wNames <- c( "wFood1", "wFood2", "wFood3", "wFood4" ) > > # AIDS > estResultA <- aidsEst( pNames, wNames, "xFood", + data = Blanciforti86[ 2:nrow( Blanciforti86 ), ], maxiter = 1, + elaFormula = "AIDS", rcovformula=1, tol=1e-7, + method = "MK:L" ) > print( diag( estResultA$ela$marshall ) ) [1] -0.4354927 -0.2350566 -0.7288220 -0.3935498 > > # LA-AIDS + formula of AIDS > estResultL1 <- aidsEst( pNames, wNames, "xFood", + data = Blanciforti86, maxiter = 1, elaFormula = "AIDS", + rcovformula=1, tol=1e-7, method = "LA:SL" ) > print( diag( estResultL1$ela$marshall ) ) [1] -0.3894355 -0.2296041 -0.7355415 -0.3272023 > > # LA-AIDS + formula of Eales + Unnevehr > estResultL2 <- aidsEst( pNames, wNames, "xFood", + data = Blanciforti86, maxiter = 1, elaFormula = "EU", + rcovformula=1, tol=1e-7, method = "LA:SL" ) > print( diag( estResultL2$ela$marshall ) ) [1] -0.6532279 -0.2002519 -0.8789371 -1.0813375 > > # LA-AIDS + formula of Chalfant: > estResultL3 <- aidsEst( pNames, wNames, "xFood", + data = Blanciforti86, maxiter = 1, elaFormula = "Ch", + rcovformula=1, tol=1e-7, method = "LA:SL" ) > print( diag( estResultL3$ela$marshall ) ) [1] -0.9808054 -0.2515941 -0.8050901 -0.7762648 > > > > cleanEx(); ..nameEx <- "aidsPx" > > ### * aidsPx > > flush(stderr()); flush(stdout()) > > ### Name: aidsPx > ### Title: Price Index for the AIDS > ### Aliases: aidsPx > ### Keywords: models > > ### ** Examples > > data( Blanciforti86 ) > aidsPx( "S", c( "pFood1", "pFood2", "pFood3", "pFood4" ), + c( "wFood1", "wFood2", "wFood3", "wFood4" ), Blanciforti86 ) [1] 4.041210 4.125092 4.084875 4.102826 4.209488 4.225516 4.218751 4.218204 [9] 4.200284 4.205226 4.237010 4.272815 4.254950 4.266069 4.281623 4.288670 [17] 4.303236 4.318542 4.337664 4.387997 4.397315 4.431085 4.479721 4.535299 [25] 4.565668 4.605170 4.743032 4.876663 4.962229 4.999500 5.063662 5.153325 > > > > cleanEx(); ..nameEx <- "aidsTestConsist" > > ### * aidsTestConsist > > flush(stderr()); flush(stdout()) > > ### Name: aidsTestConsist > ### Title: Consistency Test of the AIDS > ### Aliases: aidsTestConsist > ### Keywords: models > > ### ** Examples > > data( Blanciforti86 ) > pNames <- c( "pFood1", "pFood2", "pFood3", "pFood4" ) > wNames <- c( "wFood1", "wFood2", "wFood3", "wFood4" ) > > estResult <- aidsEst( pNames, wNames, "xFood", + data = Blanciforti86, method = "MK:L" ) > tc <- aidsTestConsist( pNames, wNames, "xFood", Blanciforti86, + coef = estResult$coef ) > tc$mPercent # great! [1] 100 > tc$cPercent # Oh, that is bad! [1] 0 > > > > cleanEx(); ..nameEx <- "heckit" > > ### * heckit > > flush(stderr()); flush(stdout()) > > ### Name: heckit > ### Title: 2-step Heckman (heckit) estimation > ### Aliases: heckit > ### Keywords: models > > ### ** Examples > > ## Greene( 2003 ): example 22.8, page 786 > data( Mroz87 ) > Mroz87$kids <- ( Mroz87$kids5 + Mroz87$kids618 > 0 ) > greene <- heckit( wage ~ exper + I( exper^2 ) + educ + city, + lfp ~ age + I( age^2 ) + faminc + kids + educ, Mroz87 ) > summary( greene ) # print summary Estimate Std. Error t value Pr(>|t|) (Intercept) -0.971247 2.059352 -0.471628 0.637436 exper 0.021061 0.062465 0.337171 0.736155 I(exper^2) 0.000137 0.001878 0.072977 0.941859 educ 0.417019 0.100255 4.159588 3.9e-05 *** city 0.443838 0.315908 1.404961 0.160769 probitLambda -1.097587 1.265848 -0.867077 0.386393 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Multiple R-Squared: 0.126368 Adjusted R-Squared: 0.116016 > summary( greene$probit ) # summary of the 1st step probit estimation Call: glm(formula = probitformula, family = binomial(link = "probit"), data = data) Deviance Residuals: Min 1Q Median 3Q Max -1.9205 -1.2261 0.7877 1.0634 1.6515 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -4.157e+00 1.404e+00 -2.961 0.003070 ** age 1.854e-01 6.621e-02 2.800 0.005107 ** I(age^2) -2.426e-03 7.762e-04 -3.125 0.001775 ** faminc 4.580e-06 4.306e-06 1.064 0.287417 kidsTRUE -4.490e-01 1.300e-01 -3.453 0.000554 *** educ 9.818e-02 2.289e-02 4.289 1.80e-05 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 1029.7 on 752 degrees of freedom Residual deviance: 981.7 on 747 degrees of freedom AIC: 993.7 Number of Fisher Scoring iterations: 4 > # this is Example 21.4, p. 681f > greene$sigma # estimated sigma [1] 3.200058 > greene$rho # estimated rho probitLambda -0.3429898 > > ## Wooldridge( 2003 ): example 17.5, page 590 > data( Mroz87 ) > wooldridge <- heckit( log( wage ) ~ educ + exper + I( exper^2 ), + lfp ~ nwifeinc + educ + exper + I( exper^2 ) + age + kids5 + kids618, Mroz87 ) > summary( wooldridge ) # summary of the 1st step probit estimation Estimate Std. Error t value Pr(>|t|) (Intercept) -0.578102 0.305007 -1.895375 0.058724 . educ 0.109065 0.015523 7.026059 0 *** exper 0.043887 0.016261 2.698911 0.007235 ** I(exper^2) -0.000859 0.000439 -1.957349 0.050964 . probitLambda 0.032261 0.133626 0.241432 0.809338 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Multiple R-Squared: 0.156935 Adjusted R-Squared: 0.148963 > # (Example 17.1, p. 562f) and 2nd step OLS regression > wooldridge$sigma # estimated sigma [1] 0.6636287 > wooldridge$rho # estimated rho probitLambda 0.04861365 > > ## example using random numbers > nObs <- 1000 > myData <- data.frame( no = c( 1:nObs ), x1 = rnorm( nObs ), x2 = rnorm( nObs ) ) > myData$y <- 2 + myData$x1 + 0.9 * rnorm( nObs ) > myData$s <- ( 2 * myData$x1 + myData$x2 + 4 * rnorm( nObs ) - 0.2 ) > 0 > myData$y[ !myData$s ] <- NA > myHeckit <- heckit( y ~ x1, s ~ x1 + x2, myData, print.level = 1 ) Estimating 1st step Probit model . . . OK Estimating 2nd step OLS model . . . OK Calculating coefficient covariance matrix . . . OK > > ## example using random numbers with IV/2SLS estimation > nObs <- 1000 > myData <- data.frame( no = c( 1:nObs ), x1 = rnorm( nObs ), x2 = rnorm( nObs ), + u = 0.5 * rnorm( nObs ) ) > myData$w <- 1 + myData$x1 + 0.2 * myData$u + 0.1 * rnorm( nObs ) > myData$y <- 2 + myData$w + myData$u > myData$s <- ( 2 * myData$x1 + myData$x2 + 4 * rnorm( nObs ) - 0.2 ) > 0 > myData$y[ !myData$s ] <- NA > myHeckit <- heckit( y ~ w, s ~ x1 + x2, data = myData, + inst = ~ x1, print.level = 1 ) Estimating 1st step Probit model . . . OK Estimating 2nd step 2SLS/IV model . . . OK Calculating coefficient covariance matrix . . . OK > > > > cleanEx(); ..nameEx <- "insertCol" > > ### * insertCol > > flush(stderr()); flush(stdout()) > > ### Name: insertCol > ### Title: Insert Column into a Matrix > ### Aliases: insertCol > ### Keywords: array > > ### ** Examples > > m <- matrix( 1:4, 2 ) > insertCol( m, 2, 5:6 ) [,1] [,2] [,3] [1,] 1 5 3 [2,] 2 6 4 > > > > cleanEx(); ..nameEx <- "insertRow" > > ### * insertRow > > flush(stderr()); flush(stdout()) > > ### Name: insertRow > ### Title: Insert Row into a Matrix > ### Aliases: insertRow > ### Keywords: array > > ### ** Examples > > m <- matrix( 1:4, 2 ) > insertRow( m, 2, 5:6 ) [,1] [,2] [1,] 1 3 [2,] 5 6 [3,] 2 4 > > > > cleanEx(); ..nameEx <- "priceIndex" > > ### * priceIndex > > flush(stderr()); flush(stdout()) > > ### Name: priceIndex > ### Title: Calculate Price Indices > ### Aliases: priceIndex > ### Keywords: models > > ### ** Examples > > data( Missong77 ) > # Laspeyres Price Indices > priceIndex( c( "p.beef", "p.veal", "p.pork" ), + c( "q.beef", "q.veal", "q.pork" ), 1, Missong77 ) 1986 1987 1988 1989 1.0000000 0.9839897 0.9966514 1.1254719 > # Paasche Price Indices > priceIndex( c( "p.beef", "p.veal", "p.pork" ), + c( "q.beef", "q.veal", "q.pork" ), 1, Missong77, "Paasche" ) 1986 1987 1988 1989 1.0000000 0.9839652 0.9966153 1.1251109 > > data( Bleymueller251 ) > # Laspeyres Price Indices > priceIndex( c( "p.A", "p.B", "p.C", "p.D" ), c("q.A", "q.B", "q.C", "q.D" ), + 1, Bleymueller251 ) 1970 1974 1978 1.000000 1.697436 1.887179 > # Paasche Price Indices > priceIndex( c( "p.A", "p.B", "p.C", "p.D" ), c("q.A", "q.B", "q.C", "q.D" ), + 1, Bleymueller251, "Paasche" ) 1970 1974 1978 1.000000 1.582278 1.799337 > > > > cleanEx(); ..nameEx <- "print.snqProfitEst" > > ### * print.snqProfitEst > > flush(stderr()); flush(stdout()) > > ### Name: print.snqProfitEst > ### Title: Print output of estimated SNQ profit function > ### Aliases: print.snqProfitEst > ### Keywords: models > > ### ** Examples > > ## Not run: library( systemfit ) > data( germanFarms ) > germanFarms$qOutput <- germanFarms$vOutput / germanFarms$pOutput > germanFarms$qVarInput <- -germanFarms$vVarInput / germanFarms$pVarInput > germanFarms$qLabor <- -germanFarms$qLabor > pNames <- c( "pOutput", "pVarInput", "pLabor" ) > qNames <- c( "qOutput", "qVarInput", "qLabor" ) > > estResult <- snqProfitEst( pNames, qNames, "land", data = germanFarms ) > print( estResult ) Estimation results of the SNQ profit function: Functional form: default (0) Estimation method: SUR Estimated coefficients: value std.err t-value prob alpha 1 -97327.4027 37623.9488 -2.5868471 1.264236e-02 alpha 2 -86797.9576 17777.5639 -4.8824439 1.114097e-05 alpha 3 -19345.1791 3775.3417 -5.1240869 4.835754e-06 beta 1 1 -38029.0727 17199.2199 -2.2110929 3.163093e-02 beta 1 2 32859.3598 15743.5034 2.0871695 4.199050e-02 beta 1 3 5169.7129 2414.6340 2.1409923 3.717372e-02 beta 2 1 32859.3598 15743.5034 2.0871695 4.199050e-02 beta 2 2 -27854.1492 14709.6934 -1.8935914 6.407267e-02 beta 2 3 -5005.2106 2293.5586 -2.1822902 3.381336e-02 beta 3 1 5169.7129 2414.6340 2.1409923 3.717372e-02 beta 3 2 -5005.2106 2293.5586 -2.1822902 3.381336e-02 beta 3 3 -164.5024 607.3175 -0.2708672 7.876083e-01 delta 1 1 10459.6184 2429.9630 4.3044352 7.790735e-05 delta 2 1 2621.1740 1081.7413 2.4231061 1.905133e-02 delta 3 1 532.5423 239.3199 2.2252319 3.060584e-02 gamma 1 1 -294.1384 115.2189 -2.5528665 1.378313e-02 R-squared values of the netput equations: qOutput qVarInput qLabor 0.9746943 0.4488089 0.3508616 Price elasticities of the netputs: pOutput pVarInput pLabor qOutput -0.3622530 0.2845257 0.077727272 qVarInput -0.5563251 0.4065936 0.149731512 qLabor -0.3365767 0.3316017 0.004974981 warning: this profit function is not convex in netput prices. > > > > cleanEx(); ..nameEx <- "quadFuncCalc" > > ### * quadFuncCalc > > flush(stderr()); flush(stdout()) > > ### Name: quadFuncCalc > ### Title: Calculate dependent variable of a quadratic function > ### Aliases: quadFuncCalc > ### Keywords: models > > ### ** Examples > > data( germanFarms ) > # output quantity: > germanFarms$qOutput <- germanFarms$vOutput / germanFarms$pOutput > # quantity of variable inputs > germanFarms$qVarInput <- germanFarms$vVarInput / germanFarms$pVarInput > # a time trend to account for technical progress: > germanFarms$time <- c(1:20) > > # estimate a quadratic production function > estResult <- quadFuncEst( "qOutput", c( "qLabor", "land", "qVarInput", "time" ), + germanFarms ) > > quadFuncCalc( c( "qLabor", "land", "qVarInput", "time" ), germanFarms, + estResult$allCoef ) 1 2 3 4 5 6 7 8 938.7012 1002.0952 1047.9269 1139.3294 1178.9317 1207.4611 1239.9975 1292.8022 9 10 11 12 13 14 15 16 1194.7349 1286.2359 1363.4601 1537.7959 1581.2848 1695.0962 1636.7475 1618.9742 17 18 19 20 1675.6702 1741.8777 1890.4868 1991.6342 > #equal to estResult$fitted > > > > cleanEx(); ..nameEx <- "quadFuncDeriv" > > ### * quadFuncDeriv > > flush(stderr()); flush(stdout()) > > ### Name: quadFuncDeriv > ### Title: Derivatives of a quadratic function > ### Aliases: quadFuncDeriv > ### Keywords: models > > ### ** Examples > > data( germanFarms ) > # output quantity: > germanFarms$qOutput <- germanFarms$vOutput / germanFarms$pOutput > # quantity of variable inputs > germanFarms$qVarInput <- germanFarms$vVarInput / germanFarms$pVarInput > # a time trend to account for technical progress: > germanFarms$time <- c(1:20) > > # estimate a quadratic production function > estResult <- quadFuncEst( "qOutput", c( "qLabor", "land", "qVarInput", "time" ), + germanFarms ) > > # compute the marginal products of the inputs > margProducts <- quadFuncDeriv( c( "qLabor", "land", "qVarInput", "time" ), + germanFarms, estResult$allCoef, estResult$allCoefCov ) > margProducts$deriv qLabor land qVarInput time 1 2546.6963 70.37321 -0.87697378 31.415311 2 2218.5034 55.78284 -1.11310030 46.483105 3 2301.6866 55.70200 -0.74397019 38.124492 4 1509.4055 52.72320 -0.89393206 43.800087 5 1398.3578 60.23202 -0.33770187 27.603970 6 1397.1054 57.99605 0.02997629 22.457520 7 1560.8988 48.97724 0.26225368 24.220550 8 1108.3537 73.05437 0.73659135 -7.177405 9 926.7102 82.48118 1.64839516 -26.829001 10 382.7553 92.52914 1.49605442 -38.131223 11 230.2419 81.91353 1.47746659 -30.332660 12 -1402.5828 82.58996 -0.62207350 -10.604369 13 -1571.1647 69.48917 -0.94239273 1.786045 14 -1089.4578 63.42216 -0.64667010 -3.018705 15 -847.2839 89.24736 0.69181322 -47.379138 16 -349.3805 89.59520 1.42135451 -61.938159 17 -106.1504 95.91727 1.40790948 -76.920291 18 115.7999 93.42042 1.34735027 -81.477061 19 366.4300 77.01233 0.79468293 -67.890172 20 525.3894 57.94714 0.22976523 -50.399087 > > > > cleanEx(); ..nameEx <- "quadFuncEst" > > ### * quadFuncEst > > flush(stderr()); flush(stdout()) > > ### Name: quadFuncEst > ### Title: Estimate a quadratic function > ### Aliases: quadFuncEst > ### Keywords: models > > ### ** Examples > > data( germanFarms ) > # output quantity: > germanFarms$qOutput <- germanFarms$vOutput / germanFarms$pOutput > # quantity of variable inputs > germanFarms$qVarInput <- germanFarms$vVarInput / germanFarms$pVarInput > # a time trend to account for technical progress: > germanFarms$time <- c(1:20) > > # estimate a quadratic production function > estResult <- quadFuncEst( "qOutput", c( "qLabor", "land", "qVarInput", "time" ), + germanFarms ) > > estResult$coef $alpha0 (Intercept) -18903.29 $alpha x1 x2 x3 x4 16390.13047 28.86098 19.22564 -150.43466 $beta [,1] [,2] [,3] [,4] [1,] -11448.03555 244.9444987 -5.29227031 -144.5663417 [2,] 244.94450 -4.7880503 -0.36717605 5.1833120 [3,] -5.29227 -0.3671761 -0.00789875 0.5062034 [4,] -144.56634 5.1833120 0.50620344 -12.9900767 > estResult$r2 [1] 0.9973797 > > > > cleanEx(); ..nameEx <- "quantityIndex" > > ### * quantityIndex > > flush(stderr()); flush(stdout()) > > ### Name: quantityIndex > ### Title: Calculate Quantity Indices > ### Aliases: quantityIndex > ### Keywords: models > > ### ** Examples > > data( Missong77 ) > # Laspeyres Quantity Indices > quantityIndex( c( "p.beef", "p.veal", "p.pork" ), + c( "q.beef", "q.veal", "q.pork" ), 1, Missong77 ) 1986 1987 1988 1989 1.0000000 0.9706006 0.8833432 0.7429413 > # Paasche Quantity Indices > quantityIndex( c( "p.beef", "p.veal", "p.pork" ), + c( "q.beef", "q.veal", "q.pork" ), 1, Missong77, "Paasche" ) 1986 1987 1988 1989 1.0000000 0.9705765 0.8833112 0.7427030 > > data( Bleymueller251 ) > # Laspeyres Quantity Indices > quantityIndex( c( "p.A", "p.B", "p.C", "p.D" ), c("q.A", "q.B", "q.C", "q.D" ), + 1, Bleymueller251 ) 1970 1974 1978 1.000000 1.215385 1.546154 > # Paasche Quantity Indices > quantityIndex( c( "p.A", "p.B", "p.C", "p.D" ), c("q.A", "q.B", "q.C", "q.D" ), + 1, Bleymueller251, "Paasche" ) 1970 1974 1978 1.000000 1.132931 1.474185 > > > > cleanEx(); ..nameEx <- "quasiconcavity" > > ### * quasiconcavity > > flush(stderr()); flush(stdout()) > > ### Name: quasiconcavity > ### Title: Test for quasiconcavity / quasiconvexity > ### Aliases: quasiconcavity quasiconvexity > ### Keywords: array > > ### ** Examples > > quasiconcavity( matrix( 0, 3, 3 ) ) [1] TRUE > > quasiconvexity( matrix( 0, 3, 3 ) ) [1] TRUE > > m <- list() > m[[1]] <- matrix( c( 0,-1,-1, -1,-2,3, -1,3,5 ), 3, 3 ) > m[[2]] <- matrix( c( 0,1,-1, 1,-2,3, -1,3,5 ), 3, 3 ) > > quasiconcavity( m ) [1] TRUE FALSE > > quasiconvexity( m ) [1] FALSE TRUE > > > > cleanEx(); ..nameEx <- "rSquared" > > ### * rSquared > > flush(stderr()); flush(stdout()) > > ### Name: rSquared > ### Title: Calculate R squared value > ### Aliases: rSquared > ### Keywords: univar multivariate array > > ### ** Examples > > data( Blanciforti86 ) > reg <- lm( wFood1 ~ pFood1 + xFood, Blanciforti86 ) > rSquared( Blanciforti86$wFood1, reg$residuals ) [,1] [1,] 0.6573009 > summary( reg )$r.squared # returns the same value [1] 0.6573009 > > > > cleanEx(); ..nameEx <- "readFront41out" > > ### * readFront41out > > flush(stderr()); flush(stdout()) > > ### Name: readFront41out > ### Title: Read output of Frontier 4.1 > ### Aliases: readFront41out > ### Keywords: models > > ### ** Examples > > data( Coelli ) > Coelli$logOutput <- log( Coelli$output ) > Coelli$logCapital <- log( Coelli$capital ) > Coelli$logLabour <- log( Coelli$labour ) > > writeFront41in( Coelli, "firm", "time", "logOutput", + c( "logCapital", "logLabour" ) ) > > ## Not run: > ##D system( "wine front41.exe" ) > ##D sfa <- readFront41out() > ##D sfa$mleResults > ##D sfa$efficiency > ##D > ## End(Not run) > > > > cleanEx(); ..nameEx <- "semidefiniteness" > > ### * semidefiniteness > > flush(stderr()); flush(stdout()) > > ### Name: semidefiniteness > ### Title: Test for negative and positive semidefiniteness > ### Aliases: semidefiniteness > ### Keywords: array > > ### ** Examples > > # a positive semidefinite matrix > semidefiniteness( matrix( 1, 3, 3 )) $negative [1] FALSE $positive [1] TRUE > > # a negative semidefinite matrix > semidefiniteness( matrix(-1, 3, 3 )) $negative [1] TRUE $positive [1] FALSE > > # a matrix that is positive and negative semidefinite > semidefiniteness( matrix( 0, 3, 3 )) $negative [1] TRUE $positive [1] TRUE > > # a matrix that is neither positive nor negative semidefinite > semidefiniteness( matrix( 1:9, 3, 3 )) $negative [1] FALSE $positive [1] FALSE > > > > cleanEx(); ..nameEx <- "snqProfitCalc" > > ### * snqProfitCalc > > flush(stderr()); flush(stdout()) > > ### Name: snqProfitCalc > ### Title: Calculations with the SNQ Profit function > ### Aliases: snqProfitCalc > ### Keywords: models > > ### ** Examples > > data( germanFarms ) > germanFarms$qOutput <- germanFarms$vOutput / germanFarms$pOutput > germanFarms$qVarInput <- -germanFarms$vVarInput / germanFarms$pVarInput > germanFarms$qLabor <- -germanFarms$qLabor > germanFarms$time <- c( 0:19 ) > pNames <- c( "pOutput", "pVarInput", "pLabor" ) > qNames <- c( "qOutput", "qVarInput", "qLabor" ) > fNames <- c( "land", "time" ) > > estResult <- snqProfitEst( pNames, qNames, fNames, data = germanFarms ) > snqProfitCalc( pNames, fNames, estResult$estData, estResult$weights, + estResult$coef ) X1 X2 X3 profit 1 93936.43 -52884.29 -12570.64 28481.49 2 98124.87 -51961.22 -12226.80 26129.94 3 102204.88 -51226.00 -11946.76 23823.26 4 107439.17 -51539.49 -12019.45 27886.48 5 111264.95 -51717.02 -12036.29 31604.50 6 114818.48 -51932.08 -12048.46 31723.63 7 120083.76 -52922.36 -12221.37 38198.87 8 122276.34 -53530.98 -12320.89 36443.14 9 120544.99 -55089.09 -12677.21 32011.30 10 124915.71 -55718.58 -12747.91 29574.55 11 128902.89 -56848.18 -12954.30 28137.87 12 144272.77 -57019.78 -12849.23 38278.31 13 150901.64 -57957.64 -12986.58 40253.68 14 157855.48 -59410.40 -13279.27 48620.90 15 156368.67 -61333.62 -13701.06 47263.18 16 156337.07 -61760.22 -13714.96 31825.31 17 165323.51 -62423.00 -13784.43 38228.85 18 171695.97 -62432.67 -13707.91 30257.78 19 182075.71 -62671.20 -13651.16 30911.46 20 192423.32 -63773.05 -13799.33 38436.82 > > > > cleanEx(); ..nameEx <- "snqProfitEla" > > ### * snqProfitEla > > flush(stderr()); flush(stdout()) > > ### Name: snqProfitEla > ### Title: Elasticities of SNQ Profit function > ### Aliases: snqProfitEla > ### Keywords: models > > ### ** Examples > > # just a stupid simple example > snqProfitEla( matrix(101:109,3,3), c(1,1,1), c(1,-1,-1), c(0.4,0.3,0.3) ) p1 p2 p3 q1 5.0 -2.20 -2.80 q2 2.4 -1.05 -1.35 q3 2.6 -1.15 -1.45 > > # now with real data > data( germanFarms ) > germanFarms$qOutput <- germanFarms$vOutput / germanFarms$pOutput > germanFarms$qVarInput <- -germanFarms$vVarInput / germanFarms$pVarInput > germanFarms$qLabor <- -germanFarms$qLabor > germanFarms$time <- c( 0:19 ) > pNames <- c( "pOutput", "pVarInput", "pLabor" ) > qNames <- c( "qOutput", "qVarInput", "qLabor" ) > > estResult <- snqProfitEst( pNames, qNames, c("land","time"), data=germanFarms ) > > estResult$ela # price elasticities at mean prices and mean quantities pOutput pVarInput pLabor qOutput 0.1704468 -0.1157277 -0.05471903 qVarInput 0.2262792 -0.1596354 -0.06664374 qLabor 0.2369458 -0.1475920 -0.08935378 > > # price elasticities at the last observation (1994/95) > snqProfitEla( estResult$coef$beta, estResult$estData[ 20, pNames ], + estResult$estData[ 20, qNames ], estResult$weights ) pOutput pVarInput pLabor qOutput 0.1460128 -0.08748643 -0.05852635 qVarInput 0.1778510 -0.11509188 -0.06275908 qLabor 0.2051378 -0.10820699 -0.09693081 > > > > cleanEx(); ..nameEx <- "snqProfitEst" > > ### * snqProfitEst > > flush(stderr()); flush(stdout()) > > ### Name: snqProfitEst > ### Title: Estimation of a SNQ Profit function > ### Aliases: snqProfitEst > ### Keywords: models > > ### ** Examples > > data( germanFarms ) > germanFarms$qOutput <- germanFarms$vOutput / germanFarms$pOutput > germanFarms$qVarInput <- -germanFarms$vVarInput / germanFarms$pVarInput > germanFarms$qLabor <- -germanFarms$qLabor > pNames <- c( "pOutput", "pVarInput", "pLabor" ) > qNames <- c( "qOutput", "qVarInput", "qLabor" ) > > estResult <- snqProfitEst( pNames, qNames, "land", data = germanFarms ) > estResult$ela # Oh, that looks bad! pOutput pVarInput pLabor qOutput -0.3622530 0.2845257 0.077727272 qVarInput -0.5563251 0.4065936 0.149731512 qLabor -0.3365767 0.3316017 0.004974981 > > # it it reasonable to account for technological progress > germanFarms$time <- c( 0:19 ) > estResult2 <- snqProfitEst( pNames, qNames, c("land","time"), data=germanFarms ) > estResult2$ela # Ah, that looks good! pOutput pVarInput pLabor qOutput 0.1704468 -0.1157277 -0.05471903 qVarInput 0.2262792 -0.1596354 -0.06664374 qLabor 0.2369458 -0.1475920 -0.08935378 > > > > cleanEx(); ..nameEx <- "snqProfitHessian" > > ### * snqProfitHessian > > flush(stderr()); flush(stdout()) > > ### Name: snqProfitHessian > ### Title: SNQ Profit function: Hessian matrix > ### Aliases: snqProfitHessian > ### Keywords: models > > ### ** Examples > > # just a stupid simple example > snqProfitHessian( matrix(101:109,3,3), c(1,1,1), c(0.4,0.3,0.3) ) [,1] [,2] [,3] [1,] 5.0 -2.20 -2.80 [2,] -2.4 1.05 1.35 [3,] -2.6 1.15 1.45 > > # now with real data > data( germanFarms ) > germanFarms$qOutput <- germanFarms$vOutput / germanFarms$pOutput > germanFarms$qVarInput <- -germanFarms$vVarInput / germanFarms$pVarInput > germanFarms$qLabor <- -germanFarms$qLabor > germanFarms$time <- c( 0:19 ) > pNames <- c( "pOutput", "pVarInput", "pLabor" ) > qNames <- c( "qOutput", "qVarInput", "qLabor" ) > > estResult <- snqProfitEst( pNames, qNames, c("land","time"), data=germanFarms ) > > estResult$hessian # the Hessian at mean prices and mean quantities pOutput pVarInput pLabor qOutput 22813.736 -12620.371 -2997.4434 qVarInput -12620.371 7254.105 1521.2194 qLabor -2997.443 1521.219 462.6152 > > # Hessian at the last observation (1994/95) > snqProfitHessian( estResult$coef$beta, estResult$estData[ 20, pNames ], + estResult$weights ) pOutput pVarInput pLabor qOutput 30240.134 -12205.214 -2973.3857 qVarInput -12205.214 5320.418 1056.5079 qLabor -2973.386 1056.508 344.6469 > > > > cleanEx(); ..nameEx <- "snqProfitHessianDeriv" > > ### * snqProfitHessianDeriv > > flush(stderr()); flush(stdout()) > > ### Name: snqProfitHessianDeriv > ### Title: SNQ Profit function: Derivatives of the Hessian > ### Aliases: snqProfitHessianDeriv > ### Keywords: models > > ### ** Examples > > # just a stupid simple example > snqProfitHessianDeriv( c(1,2,3),c(0.4,0.3,0.3) ) [,1] [,2] [,3] [,4] [,5] [,6] [1,] 0 0 0 1.06283715 0.3149147 0.02332702 [2,] 0 0 0 0.23618603 0.9010060 0.12829859 [3,] 0 0 0 0.05248579 0.3848958 0.70564222 > > # now with real data > data( germanFarms ) > germanFarms$qOutput <- germanFarms$vOutput / germanFarms$pOutput > germanFarms$qVarInput <- -germanFarms$vVarInput / germanFarms$pVarInput > germanFarms$qLabor <- -germanFarms$qLabor > germanFarms$time <- c( 0:19 ) > pNames <- c( "pOutput", "pVarInput", "pLabor" ) > qNames <- c( "qOutput", "qVarInput", "qLabor" ) > > estResult <- snqProfitEst( pNames, qNames, c("land","time"), data=germanFarms ) > > snqProfitHessianDeriv( estResult$pMean, estResult$weights, 2 ) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [1,] 0 0 0 2.8517975 2.1771146 0.4155123 0 0 0 0 0 [2,] 0 0 0 0.5384073 2.2188259 0.7684992 0 0 0 0 0 [3,] 0 0 0 0.1016490 0.7602092 1.4213564 0 0 0 0 0 [,12] [,13] [,14] [,15] [1,] 0 0 0 0 [2,] 0 0 0 0 [3,] 0 0 0 0 > > > > cleanEx(); ..nameEx <- "snqProfitImposeConvexity" > > ### * snqProfitImposeConvexity > > flush(stderr()); flush(stdout()) > > ### Name: snqProfitImposeConvexity > ### Title: Imposing Convexity on a SNQ Profit function > ### Aliases: snqProfitImposeConvexity > ### Keywords: models > > ### ** Examples > > data( germanFarms ) > germanFarms$qOutput <- germanFarms$vOutput / germanFarms$pOutput > germanFarms$qVarInput <- -germanFarms$vVarInput / germanFarms$pVarInput > germanFarms$qLabor <- -germanFarms$qLabor > germanFarms$time <- c( 0:19 ) > pNames <- c( "pOutput", "pVarInput", "pLabor" ) > qNames <- c( "qOutput", "qVarInput", "qLabor" ) > fNames <- c( "land", "time" ) > estResult <- snqProfitEst( pNames, qNames, fNames, data = germanFarms ) > estResult # Note: it is NOT convex in netput prices Estimation results of the SNQ profit function: Functional form: default (0) Estimation method: SUR Estimated coefficients: value std.err t-value prob alpha 1 -79970.1603 115640.3820 -0.6915418 0.49277769 alpha 2 -85851.7559 50774.6633 -1.6908385 0.09778375 alpha 3 -20576.0170 11157.9147 -1.8440737 0.07176223 beta 1 1 16779.6820 24973.7900 0.6718917 0.50508772 beta 1 2 -13637.9458 23692.2568 -0.5756288 0.56773369 beta 1 3 -3141.7362 3092.5983 -1.0158889 0.31511342 beta 2 1 -13637.9458 23692.2568 -0.5756288 0.56773369 beta 2 2 11197.8664 22868.5218 0.4896629 0.62675112 beta 2 3 2440.0794 2881.5153 0.8468043 0.40158531 beta 3 1 -3141.7362 3092.5983 -1.0158889 0.31511342 beta 3 2 2440.0794 2881.5153 0.8468043 0.40158531 beta 3 3 701.6568 637.1537 1.1012362 0.27664953 delta 1 1 13791.9194 12077.5455 1.1419472 0.25951605 delta 1 2 -8643.8822 11738.0893 -0.7363960 0.46531148 delta 2 1 4059.0621 5083.0146 0.7985541 0.42874354 delta 2 2 -5696.8272 5113.4235 -1.1140926 0.27115475 delta 3 1 943.2179 1145.5429 0.8233807 0.41463398 delta 3 2 -1261.9766 1133.4696 -1.1133749 0.27145944 gamma 1 1 -821.4802 937.1857 -0.8765394 0.38539227 gamma 1 2 923.8323 860.8056 1.0732183 0.28889477 gamma 2 1 923.8323 860.8056 1.0732183 0.28889477 gamma 2 2 -1063.7748 684.6548 -1.5537388 0.12725228 R-squared values of the netput equations: qOutput qVarInput qLabor 0.9880570 0.5018285 0.6090835 Price elasticities of the netputs: pOutput pVarInput pLabor qOutput 0.1704468 -0.1157277 -0.05471903 qVarInput 0.2262792 -0.1596354 -0.06664374 qLabor 0.2369458 -0.1475920 -0.08935378 > estResultConvex <- snqProfitImposeConvexity( estResult ) > estResultConvex # now it is convex Estimation results of the SNQ profit function with convexityimposed: Functional form: default (0) Minimum distance estimation: convergence achieved after 224 iterations Estimated coefficients: [,1] alpha 1 -79970.184 alpha 2 -85852.379 alpha 3 -20576.020 beta 1 1 16782.172 beta 1 2 -13640.526 beta 1 3 -3141.646 beta 2 1 -13640.526 beta 2 2 11200.558 beta 2 3 2439.968 beta 3 1 -3141.646 beta 3 2 2439.968 beta 3 3 701.678 delta 1 1 13792.034 delta 1 2 -8644.080 delta 2 1 4059.137 delta 2 2 -5696.985 delta 3 1 943.228 delta 3 2 -1261.999 gamma 1 1 -821.495 gamma 1 2 923.852 gamma 2 1 923.852 gamma 2 2 -1063.792 R-squared values of the netput equations: qOutput qVarInput qLabor 0.988057 0.501831 0.609083 Price elasticities of the netputs: pOutput pVarInput pLabor qOutput 0.170466 -0.115749 -0.0547172 qVarInput 0.226321 -0.159682 -0.0666387 qLabor 0.236938 -0.147581 -0.0893573 > > > > cleanEx(); ..nameEx <- "snqProfitShadowPrices" > > ### * snqProfitShadowPrices > > flush(stderr()); flush(stdout()) > > ### Name: snqProfitShadowPrices > ### Title: Shadow Prices of a SNQ Profit function > ### Aliases: snqProfitShadowPrices > ### Keywords: models > > ### ** Examples > > data( germanFarms ) > germanFarms$qOutput <- germanFarms$vOutput / germanFarms$pOutput > germanFarms$qVarInput <- -germanFarms$vVarInput / germanFarms$pVarInput > germanFarms$qLabor <- -germanFarms$qLabor > germanFarms$time <- c( 0:19 ) > pNames <- c( "pOutput", "pVarInput", "pLabor" ) > qNames <- c( "qOutput", "qVarInput", "qLabor" ) > fNames <- c( "land", "time" ) > > estResult <- snqProfitEst( pNames, qNames, fNames, data = germanFarms ) > > snqProfitShadowPrices( pNames, fNames, estResult ) land time 1 -428.43639 6014.9903 2 -355.35459 5981.6045 3 81.35065 5228.0054 4 868.21463 4116.4051 5 1918.70282 3098.3041 6 2660.02047 2213.5511 7 3324.35549 1846.7244 8 4388.12981 531.4822 9 6350.86338 -1704.4716 10 6463.74284 -2107.7581 11 6598.02443 -2573.2064 12 4423.58217 -468.9690 13 4245.09209 -414.7765 14 4480.66062 -469.6117 15 6482.37804 -2516.0135 16 6900.81745 -3545.9529 17 6628.74029 -3168.8668 18 5762.60384 -2800.3789 19 4523.30397 -1702.8557 20 3733.52448 -734.6891 > > > > cleanEx(); ..nameEx <- "snqProfitWeights" > > ### * snqProfitWeights > > flush(stderr()); flush(stdout()) > > ### Name: snqProfitWeights > ### Title: SNQ Profit function: Weights of prices for normalization > ### Aliases: snqProfitWeights > ### Keywords: models > > ### ** Examples > > data( germanFarms ) > germanFarms$qOutput <- germanFarms$vOutput / germanFarms$pOutput > germanFarms$qVarInput <- -germanFarms$vVarInput / germanFarms$pVarInput > germanFarms$qLabor <- -germanFarms$qLabor > pNames <- c( "pOutput", "pVarInput", "pLabor" ) > qNames <- c( "qOutput", "qVarInput", "qLabor" ) > snqProfitWeights( pNames, qNames, germanFarms ) [1] 0.66172131 0.27573702 0.06254167 > > > > cleanEx(); ..nameEx <- "summary.heckit" > > ### * summary.heckit > > flush(stderr()); flush(stdout()) > > ### Name: summary.heckit > ### Title: Print output of heckit estimation > ### Aliases: summary.heckit > ### Keywords: models > > ### ** Examples > > data( Mroz87 ) > wooldridge <- heckit( log( wage ) ~ educ + exper + I( exper^2 ), + lfp ~ nwifeinc + educ + exper + I( exper^2 ) + age + kids5 + kids618, Mroz87 ) > summary( wooldridge ) # summary of the 1st step probit estimation Estimate Std. Error t value Pr(>|t|) (Intercept) -0.578102 0.305007 -1.895375 0.058724 . educ 0.109065 0.015523 7.026059 0 *** exper 0.043887 0.016261 2.698911 0.007235 ** I(exper^2) -0.000859 0.000439 -1.957349 0.050964 . probitLambda 0.032261 0.133626 0.241432 0.809338 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Multiple R-Squared: 0.156935 Adjusted R-Squared: 0.148963 > # (Example 17.1, p. 562f) and 2nd step OLS regression > > > > cleanEx(); ..nameEx <- "translogCalc" > > ### * translogCalc > > flush(stderr()); flush(stdout()) > > ### Name: translogCalc > ### Title: Calculate dependent variable of a translog function > ### Aliases: translogCalc > ### Keywords: models > > ### ** Examples > > data( germanFarms ) > # output quantity: > germanFarms$qOutput <- germanFarms$vOutput / germanFarms$pOutput > # quantity of variable inputs > germanFarms$qVarInput <- germanFarms$vVarInput / germanFarms$pVarInput > # a time trend to account for technical progress: > germanFarms$time <- c(1:20) > > # estimate a quadratic production function > estResult <- translogEst( "qOutput", c( "qLabor", "land", "qVarInput", "time" ), + germanFarms ) > > translogCalc( c( "qLabor", "land", "qVarInput", "time" ), germanFarms, + estResult$allCoef ) 1 2 3 4 5 6 7 8 943.9667 975.9852 1060.6385 1145.4517 1188.1554 1206.2718 1236.5526 1292.7237 9 10 11 12 13 14 15 16 1196.9841 1281.0442 1362.2633 1542.6942 1575.6501 1699.6713 1632.7159 1618.6759 17 18 19 20 1678.3680 1741.8183 1887.9169 1992.8400 > #equal to estResult$fitted > > > > cleanEx(); ..nameEx <- "translogDeriv" > > ### * translogDeriv > > flush(stderr()); flush(stdout()) > > ### Name: translogDeriv > ### Title: Derivatives of a translog function > ### Aliases: translogDeriv > ### Keywords: models > > ### ** Examples > > data( germanFarms ) > # output quantity: > germanFarms$qOutput <- germanFarms$vOutput / germanFarms$pOutput > # quantity of variable inputs > germanFarms$qVarInput <- germanFarms$vVarInput / germanFarms$pVarInput > # a time trend to account for technical progress: > germanFarms$time <- c(1:20) > > # estimate a quadratic production function > estResult <- translogEst( "qOutput", c( "qLabor", "land", "qVarInput", "time" ), + germanFarms ) > > # compute the marginal products of the inputs > margProducts <- translogDeriv( c( "qLabor", "land", "qVarInput", "time" ), + germanFarms, estResult$allCoef, estResult$allCoefCov ) > margProducts$deriv qLabor land qVarInput time 1 6137.01012 163.84424 -7.0437159 208.28371 2 4381.94904 123.58899 -4.3915969 110.90773 3 4012.40673 107.02978 -2.8183425 47.71083 4 2675.93516 105.71777 -2.3880325 45.81646 5 2202.67474 98.88737 -1.2879319 18.53370 6 1859.76689 82.11185 -0.4312382 16.62160 7 1804.01073 60.26294 0.2283467 23.31485 8 1340.31183 81.26647 0.8116562 -17.58298 9 727.94966 68.14356 1.7748328 -14.91299 10 291.56187 86.04333 1.5521664 -33.34425 11 -44.49238 74.82655 1.5188881 -21.32448 12 -1280.40242 126.38542 -1.3315384 -47.81918 13 -1527.38834 114.99121 -1.6257321 -38.28692 14 -1105.00422 103.10888 -1.2105681 -39.12230 15 -741.70138 102.00035 0.2458966 -53.68207 16 -189.44046 82.85026 1.2623097 -51.03273 17 453.37590 87.66300 1.4376562 -66.31119 18 920.19258 84.04253 1.4892366 -70.92061 19 1400.46605 76.03035 0.9669266 -69.91835 20 1662.32190 64.62612 0.3283968 -63.22596 > > > > cleanEx(); ..nameEx <- "translogEst" > > ### * translogEst > > flush(stderr()); flush(stdout()) > > ### Name: translogEst > ### Title: Estimate a translog function > ### Aliases: translogEst > ### Keywords: models > > ### ** Examples > > data( germanFarms ) > # output quantity: > germanFarms$qOutput <- germanFarms$vOutput / germanFarms$pOutput > # quantity of variable inputs > germanFarms$qVarInput <- germanFarms$vVarInput / germanFarms$pVarInput > # a time trend to account for technical progress: > germanFarms$time <- c(1:20) > > # estimate a quadratic production function > estResult <- translogEst( "qOutput", c( "qLabor", "land", "qVarInput", "time" ), + germanFarms ) > > estResult$coef $alpha0 (Intercept) -133.4852 $alpha x1 x2 x3 x4 24.58924 19.92287 34.63080 -12.46742 $beta [,1] [,2] [,3] [,4] [1,] -16.864738 14.3472829 -8.591123 -3.5734031 [2,] 14.347283 0.2644969 -3.468479 -1.1633161 [3,] -8.591123 -3.4684787 -3.909638 2.7933263 [4,] -3.573403 -1.1633161 2.793326 -0.2456535 > estResult$r2 [1] 0.9982506 > > > > cleanEx(); ..nameEx <- "translogHessian" > > ### * translogHessian > > flush(stderr()); flush(stdout()) > > ### Name: translogHessian > ### Title: Hessian matrix of a translog function > ### Aliases: translogHessian > ### Keywords: models > > ### ** Examples > > data( germanFarms ) > # output quantity: > germanFarms$qOutput <- germanFarms$vOutput / germanFarms$pOutput > # quantity of variable inputs > germanFarms$qVarInput <- germanFarms$vVarInput / germanFarms$pVarInput > # a time trend to account for technical progress: > germanFarms$time <- c(1:20) > > # estimate a quadratic production function > estResult <- translogEst( "qOutput", c( "qLabor", "land", "qVarInput", "time" ), + germanFarms ) > > # compute the Hessian matrices > hessians <- translogHessian( c( "qLabor", "land", "qVarInput", "time" ), + germanFarms, estResult$allCoef ) > hessians[[ 1 ]] [,1] [,2] [,3] [,4] [1,] 23383.897716 1467.12823 -55.912414 -988.366735 [2,] 1467.128231 -3805.06084 -1.473985 -10.776912 [3,] -55.912414 -988.36673 3920.200523 3.183619 [4,] -1.473985 -10.77691 3.183619 -394.215222 > > # compute the bordered Hessian matrices > borderedHessians <- translogHessian( c( "qLabor", "land", "qVarInput", "time" ), + germanFarms, estResult$allCoef, bordered = TRUE ) > borderedHessians[[ 1 ]] [,1] [,2] [,3] [,4] [,5] [1,] 0.000000 6137.010120 163.84424 -7.043716 208.283709 [2,] 6137.010120 23383.897716 1467.12823 -55.912414 -988.366735 [3,] 163.844242 -55.912414 -3805.06084 -1.473985 -10.776912 [4,] 1467.128231 -1.473985 -988.36673 3920.200523 3.183619 [5,] -7.043716 208.283709 -10.77691 3.183619 -394.215222 > > > > cleanEx(); ..nameEx <- "triang" > > ### * triang > > flush(stderr()); flush(stdout()) > > ### Name: triang > ### Title: Upper triangular matrix from a vector > ### Aliases: triang > ### Keywords: array > > ### ** Examples > > v <- c( 1:5 ) > triang( v, 3 ) [,1] [,2] [,3] [1,] 1 2 3 [2,] 0 4 5 [3,] 0 0 0 > > > > cleanEx(); ..nameEx <- "vecli" > > ### * vecli > > flush(stderr()); flush(stdout()) > > ### Name: vecli > ### Title: Vector of linear independent values > ### Aliases: vecli > ### Keywords: array > > ### ** Examples > > # a symmetric n x n matrix > m <- cbind(c(11,12,13),c(12,22,23),c(13,23,33)) > vecli(m) # returns: 11 12 13 22 23 33 [1] 11 12 13 22 23 33 > > > > cleanEx(); ..nameEx <- "vecli2m" > > ### * vecli2m > > flush(stderr()); flush(stdout()) > > ### Name: vecli2m > ### Title: Convert vector of linear independent values into a Matrix > ### Aliases: vecli2m > ### Keywords: array > > ### ** Examples > > v <- c( 11, 12, 13, 22, 23, 33 ) > vecli2m( v ) [,1] [,2] [,3] [1,] 11 12 13 [2,] 12 22 23 [3,] 13 23 33 > > > > cleanEx(); ..nameEx <- "veclipos" > > ### * veclipos > > flush(stderr()); flush(stdout()) > > ### Name: veclipos > ### Title: Position in a vector of linear independent values > ### Aliases: veclipos > ### Keywords: array > > ### ** Examples > > veclipos( 1, 2, 3 ) # returns: 2 [1] 2 > > > > cleanEx(); ..nameEx <- "writeFront41in" > > ### * writeFront41in > > flush(stderr()); flush(stdout()) > > ### Name: writeFront41in > ### Title: Write instruction file for Frontier 4.1 > ### Aliases: writeFront41in > ### Keywords: models > > ### ** Examples > > data( Coelli ) > Coelli$logOutput <- log( Coelli$output ) > Coelli$logCapital <- log( Coelli$capital ) > Coelli$logLabour <- log( Coelli$labour ) > > writeFront41in( Coelli, "firm", "time", "logOutput", + c( "logCapital", "logLabour" ) ) > > ## Not run: > ##D system( "wine front41.exe" ) > ##D > ## End(Not run) > > > > ### *