set echo off
set messages off

function series GDV(series x, scalar m, scalar s)
     series den = exp(pi*(m-x)/(s*sqrt(3)))
     return 1/(1+den)
end function

open gretl_plant.gdt -q

# CHOICE here: do we want fixed effects or not?
scalar FE = 1
# scalar FE = 0

scalar SSRmin = 1e10
scalar SSRmax = 0
scalar optmu
scalar optsig
scalar badmu
scalar badsig
scalar nvals = 0

# apply scaling (though not for HHI)
Price /=100
Domestic_Sales /= 10000
Export_Sales /= 10000
Total_import /= 10000

# range and step for mu
scalar m_min = 0.28
scalar m_max = 0.501
scalar m_step = 0.01

# range and step for sigma
scalar s_min = 0.01
scalar s_max = 0.201
scalar s_step = 0.01

# regular regressors: include Cost?
list X = const Cost Domestic_Sales Export_Sales Soda PCFT Total_import
# logistic term
series Log
# dependent variable
series y

/**** GRID SEARCH ****/

loop for (mu=m_min; mu<m_max; mu+=m_step) --quiet
  loop for (sig=s_min; sig<s_max; sig+=s_step) --quiet
    Log = GDV(HHI, mu, sig)
    y = Price - Log
    if FE
      panel y X --quiet
    else
      ols y X --quiet
    endif
    SSR = $ess
    if SSR < SSRmin
      optmu = mu
      optsig = sig
      SSRmin = SSR
    endif
    if SSR > SSRmax
      badmu = mu
      badsig = sig
      SSRmax = SSR
    endif    
    nvals++
  endloop
endloop

printf "*** Estimated via GRID SEARCH\n\n"
printf "mu in [%g,%g), sigma in [%g,%g)\n", m_min, m_max, 
  s_min, s_max
printf "number of points evaluated: %d\n", nvals
printf "min(SSR): %g for mu = %.3f, sigma = %.3f\n",
  SSRmin, optmu, optsig
printf "max(SSR): %g for mu = %.3f, sigma = %.3f\n",
  SSRmax, badmu, badsig

Log = GDV(HHI, optmu, optsig)
y = Price - Log
if FE
  panel y X
else
  ols y X
endif
  
/**** NLS ****/

# initialize NLS from grid search results
scalar m = optmu
scalar s = optsig

printf "*** Estimated via NLS\n"

if FE
  series u = $unit
  list D = dummify(u)
  matrix b = $coeff
  matrix g = zeros(nelem(D),1)
  nls Price = lincomb(X,b) + lincomb(D,g) + Log
     series Log = GDV(HHI, m, s)
     params b g m s
  end nls -q 
  b |= g
  vnam = varname(X) ~ "," ~ varname(D) 
else
  matrix b = $coeff
  nls Price = lincomb(X,b) + Log
     series Log = GDV(HHI, m, s)
     params b m s
  end nls -q
  vnam = varname(X)
endif

# print out NLS results with nice formatting
matrix cf = b | {m; s}
matrix se = sqrt(diag($vcv))
vnam = vnam ~ ",mu,sigma"
cfse = cf ~ se
modprint cfse vnam
printf "SSR = %g, R-squared = %g\n\n", $ess, $rsq

if m < m_min || m > m_max
  printf "NLS: mu (%g) is not in [%g,%g)\n", m, m_min, m_max
endif
if s < s_min || s > s_max
  printf "NLS: sigma (%g) is not in [%g,%g)\n", s, s_min, s_max
endif