Announcement

Collapse
No announcement yet.
X
  • Filter
  • Time
  • Show
Clear All
new posts

  • Median Unbiased Estimators Laubach Williams

    Hello everyone

    I have been trying to replicate the Laubach and Williams paper: Measuring the Natural Rate of Interest. Link:http://www.frbsf.org/economic-resear.../wp2016-11.pdf
    I have already been able to compute the first part of the estimation (obtain lambda_g) but I don't know how to impose it on the second stage of the model. Details of the procedure are given in page 7.

    My question is: how do I tell Stata that there is a specific value for the ratio between the standard deviations of the state variables? Is it on the constraints? Or do I need to create a specific state equation just for this?

    Thank you in advance
    Greetings
    Alexandre Mendonça

  • #2
    If there is no obvious way of getting the code the authors used in the paper, then I would kindly ask the authors for a copy of the code so that you can reproduce the results. Good luck! If they provide you the code, then Stata the task is much easier. Once you can formulate in math (formulas, not words) what you want Stata to do, and possibly provide the SAS/R/whatever code then you can re-post here for Stata advice.

    Comment


    • #3
      The authors do provide the code but it is written in R, which I am not familiar with.
      Please find the code attached in the following link http://www.frbsf.org/economic-resear...john-williams/ under supplement (+)

      Comment


      • #4
        This is the R code for the second stage that you want to replicate in Stata:
        Code:
        #------------------------------------------------------------------------------#
        # File:        median.unbiased.estimator.stage2.R
        #
        # Description: This file implements the median unbiased estimation of the
        #              signal-to-noise ratio lambda_z following Stock and Watson (1998).
        #------------------------------------------------------------------------------#
        median.unbiased.estimator.stage2 <- function(y, x) {
            T <- dim(x)[1]
            stat <- rep(0, T-2*4+1)
            for (i in 4:(T-4)) {
                xr <- cbind(x, c(rep(0, i), rep(1, T-i)))
                xi <- solve(t(xr)%*%xr)
                b <- solve(t(xr)%*%xr,t(xr)%*%y)
                s3 <- sum((y-xr%*%b)^2)/(T-dim(xr)[2])
                stat[i+1-4] <- b[dim(xr)[2]]/sqrt(s3*xi[dim(xr)[2],dim(xr)[2]])
            }
            ew <- 0
            for (i in 1:length(stat)) {
                ew <- ew+exp((stat[i]^2)/2)
            }
            ew <- log(ew/length(stat))
            mw <- mean(stat^2)
            qlr <- max(stat^2)
        
            # Values are from Table 3 in Stock and Watson (1998)
            # Test Statistic: Exponential Wald (EW)
            valew <- c(0.426, 0.476, 0.516, 0.661, 0.826, 1.111,
                       1.419, 1.762, 2.355, 2.91,  3.413, 3.868, 4.925,
                       5.684, 6.670, 7.690, 8.477, 9.191, 10.693, 12.024,
                       13.089, 14.440, 16.191, 17.332, 18.699, 20.464,
                       21.667, 23.851, 25.538, 26.762, 27.874)
            # Test Statistic: Mean Wald (MW)
            valmw <- c(0.689, 0.757, 0.806, 1.015, 1.234, 1.632,
                       2.018, 2.390, 3.081, 3.699, 4.222, 4.776, 5.767,
                       6.586, 7.703, 8.683, 9.467, 10.101, 11.639, 13.039,
                       13.900, 15.214, 16.806, 18.330, 19.020, 20.562,
                       21.837, 24.350, 26.248, 27.089, 27.758)
            # Test Statistic: QLR
            valql <- c(3.198, 3.416, 3.594, 4.106, 4.848, 5.689,
                       6.682, 7.626, 9.16,  10.66, 11.841, 13.098, 15.451,
                       17.094, 19.423, 21.682, 23.342, 24.920, 28.174, 30.736,
                       33.313, 36.109, 39.673, 41.955, 45.056, 48.647, 50.983,
                       55.514, 59.278, 61.311, 64.016)
            
            lame <- NA
            lamm <- NA
            lamq <- NA
            
            # Median-unbiased estimator of lambda_g for given values of the test
            # statistics are obtained using the procedure described in the 
            # footnote to Stock and Watson (1998) Table 3.
            if (ew <= valew[1]) {
                lame <- 0
            } else {
                for (i in 1:(length(valew)-1)) {
                    if ((ew > valew[i]) & (ew <= valew[i+1])) {
                        lame <- i-1+(ew-valew[i])/(valew[i+1]-valew[i])
                    }
                }
            }
        
            if (mw <= valmw[1]) {
                lamm <- 0
            } else {
                for (i in 1:(length(valmw)-1)) {
                    if ((mw > valmw[i]) & (mw <= valmw[i+1])) {
                        lamm <- i-1+(mw-valmw[i])/(valmw[i+1]-valmw[i])
                    }
                }
            }
        
            if (qlr <= valql[1]) {
                lamq <- 0
            } else {
                for (i in 1:(length(valql)-1)) {
                    if ((qlr > valql[i]) & (qlr <= valql[i+1])) {
                        lamq <- i-1+(qlr-valql[i])/(valql[i+1]-valql[i])
                    }
                }
            }
            if (is.na(lame) | is.na(lamm) | is.na(lamq)) {
                    print("At least one statistic has an NA value. Check to see if your EW, MW, and/or QLR value is outside of Table 3.")
            }
            
            stats <- c(ew,mw,qlr)
            lams  <- c(lame,lamm,lamq)
            return(lame/T)
        }
        I would first replicate the analysis in R to verify the R code runs without errors, and to get the R output. It does not require much R skills to run existing code. You can install R from here: https://cran.r-project.org/. I recommend you also install RStudio. It is possible that you first have to replicate stage 1 but the code seems self-contained to me. Personally, I am not interested in doing this but I might be intefested more if you can narrow down the problem to parts of the R code.

        What is the output of the R code above, and which part do you need help with to translate into Stata?

        Comment


        • #5
          The authors also provide a code guide, "HLW_Code_Guide.pdf", in the same supplement folder. The authors say that you need to run "run.hlw.R", so please do that first. Then, please report the R output from stage 2 so that we know what the correct result R and Stata should get.

          Comment


          • #6
            By looking at the Code Guide my doubt is related with matrix Q in section 7.4 of page 10. How do I impose such a covariance matrix on stata?

            Thank you for all the help

            Comment


            • #7
              The supplemental file unpack.parameters.stage2.R gives the R code for the covariance matrix:

              Code:
                Q         <- matrix(0, 4, 4)
                Q[1, 1]   <- parameters[10]^2              # sigma_4
                Q[4, 4]   <- (lambda.g * parameters[10])^2 # sigma_4
              Stata has two matrix programming languages, the older matrix language and Mata. There is also a user-written program Rcall (GitHub) which supposedly now can construct an identical matrix in Stata, as announced here: http://www.statalist.org/forums/foru...40-rcall-1-5-0. I will only illustrate the first of these three methods. In Stata (not Mata), the basic syntax for any command is

              matrix name = matrix-expression

              So, to create a similar 4*4 matrix Q, type
              Code:
              matrix Q = (1,0,0,0\0,0,0,0\0,0,0,0\0,0,0,16)
              To list the matrix, type

              Code:
              matrix list Q
              Stata will output

              Code:
              symmetric Q[4,4]
                  c1  c2  c3  c4
              r1   1
              r2   0   0
              r3   0   0   0
              r4   0   0   0  16
              I am not sure what values you want instead of 1 and 16 because I do not know the values of "sigma_4"" and "lambda.g". Whether you actually have to create the matrix Q is another issue. I answered according to my ability. My response should give you an overview of what is involved. I am not good in matrix programming but others on Statalist are.

              Comment


              • #8
                First of all, thank you for this contribution, I really appreciate it!
                I know how to create a matrix but the problem here is that the sspace command does not allow me to set an initial covariance matrix, it only allows me to chose its structure (diagonal, identity or unstructured).
                My code is the following:

                constraint 1 [gdpnat]L.gdpnat = 1
                constraint 2 [gdpnat]L.growth = 1
                constraint 3 [gdpnat_l1]L.gdpnat = 1
                constraint 4 [gdpnat_l2]L.gdpnat_l1 = 1
                constraint 5 [growth]L.growth = 1
                constraint 6 [growth_l1]L.growth = 1
                constraint 7 [growth_l2]L.growth_l1 = 1
                constraint 8 [y]gdpnat = 1
                constraint 9 [y]gdpnat_l1 = -[y]L.y
                constraint 10 [y]gdpnat_l2 = -[y]L2.y
                constraint 11 [y]growth_l1 = -[y]L.ffr
                constraint 12 [y]growth_l2 = -[y]L2.ff
                constraint 13 [infl]L.y = -[infl]gdpnat_l1
                constraint 14 [infl]L2.infl = (1 - [infl]L.infl)/3
                constraint 15 [infl]L3.infl = [infl]L2.infl
                constraint 16 [infl]L4.infl = [infl]L2.infl


                sspace (gdpnat L.gdpnat L.growth, state noconstant) ///
                (gdpnat_l1 L.gdpnat, state noconstant noerror) ///
                (gdpnat_l2 L.gdpnat_l1, state noconstant noerror) ///
                (growth L.growth, state noconstant) ///
                (growth_l1 L.growth, state noconstant noerror) ///
                (growth_l2 L.growth_l1, state noconstant noerror) ///
                (y gdpnat L.y gdpnat_l1 L2.y gdpnat_l2 L.ffr L2.ffr growth_l1 growth_l2, ) ///
                (infl L.infl L2.infl L3.infl L4.infl L.y gdpnat_l1, noconstant), ///
                constraints(1/16) covobserved(diagonal) from(init2)

                init2 is the initial matrix for this model. It is obtained from OLS regressions identical to each of these regressions. Therefore, I am setting a value for both the variance of the gdpnat state equation and for the variance of the growth state equation. How do I impose a relationship between both of these variances?
                Please note that the matrix init2 does not include only variances and covariances, it also includes the OLS estimates of the coefficients (betas and alphas).

                Thank you in advance

                Comment


                • #9
                  Please note that the values of sigma_4 and lambda_g are not the main issue here. My question is how to impose a value for the relationship between standard errors of different state equations

                  Comment


                  • #10
                    Ok, and welcome to Statalist. Next time, please state in the original posting which Stata command you use. Use the forum [CODE] tag for showing code. Sorry, I cannot help more on this topic.

                    Comment


                    • #11
                      Hi Alexandre, were you able to resolve this issue? Could you also share your process of obtaining the first stage lambda estimate?

                      Thank you!

                      Comment

                      Working...
                      X