```
string lbl = argname(x)
if method == 0 # do nothing
series ret = x
elif method == 1 # get previous value
genr time
series OK = ok(x)
series ret = x
series tmp = OK ? time : NA
scalar ini = min(tmp)
scalar fin = max(tmp)
smpl ini fin
ret = OK ? x : ret(-1)
string lbl = sprintf("gap-filled version of %s (with repetition)", argname(x))
setinfo ret --description="@lbl"
elif method == 2 # interpolate_linearly
set warnings off
series ret = lin_int(x)
string lbl = sprintf("gap-filled version of %s (with interpolation)", argname(x))
setinfo ret --description="@lbl"
endif
return ret
```

```
# phi (as in p-high) is optional and defines an upper tail mass
# different from lower (the default phi == 0 means ignore)
#
# This is a rewrite of the function in the winsor.gfn package
# (originally by JoshuaHe2015@163.com)
smpl --no-missing x
# standard symmetric or asymmetric case?
phi = !phi ? 1 - p : phi
# determine bounds
matrix lowhi = quantile({x}, {p, phi})
# lower end (and address non-existing extreme quantiles)
scalar low = ok(lowhi[1]) ? lowhi[1] : min(x)
x = (x < low) ? low : x
# upper end
scalar hi = ok(lowhi[2]) ? lowhi[2] : max(x)
x = (x > hi) ? hi : x
# prepare output
string label = sprintf("winsorized %s (%g,%g)", argname(x), p, phi)
setinfo x --description="@label"
return x
```

```
# FIXME: cover the case with some variances == 0
# (use misszero() after scaling, or something)
#
# Forces the matrix m into the positive semi-definite region.
#
# Ported from "DomPazz" in Stackoverflow, apparently
# mimicking the nearPD() function in R.
# Because of re-scaling ( to correlation matrix ), the
# epsilon should implicitly apply to the correlation-based
# eigenvalues.
#
# The return value 0 or 1 indicates whether m was altered or not.
matrix s = sqrt(diag(m)) # std dev
matrix scaling = s * s'
matrix ms = m ./ scaling # rescaled
matrix eigvec
matrix eigval = eigensym(ms, &eigvec)
matrix val = (eigval .> epsilon) ? eigval : epsilon # like xmax
if sum(val .> eigval) # actually something was changed
matrix T = 1 / ( (eigvec .^2) * val )
# transform vector T to a diagonal matrix
matrix temp = I(rows(T))
temp[diag] = sqrt(T)
# and also val
matrix temp2 = I(rows(T))
temp2[diag] = sqrt(val)
matrix B = temp * eigvec * temp2
ms = B * B'
# undo the re-scaling
m = ms .* scaling
return 1
else
return 0
endif
```

```
# Sets elements to zero if they are really close.
# The return value 0 or 1 indicates whether m was altered or not.
# (an older version copied and returned the matrix)
matrix indicator = (abs(m) .< thresh)
m = indicator ? 0 : m
if sum(indicator)
return 1
else
return 0
endif
```

```
# Transforms comma-separated string to an array of strings.
# (The comma as separation character can be overridden with
# the sep argument, but only the first character is used.)
if !exists(sep)
string sep = ","
else
sep = substr(sep,1,1) # only first char
endif
strings out = strsplit(strsub(in, sep, " "))
return out
```

```
# returns A premultiplied by K_mn,
# K_mn being the commutation matrix, see Magnus/Neudecker
# (more efficient than explicit pre-multiplication,
# just reshuffles stuff)
# If post != 0, does post-multiplication:
# if n == 0 (default), then it's understood to be equal to m,
# as per Magnus/Neudecker convention: K_nn = K_n
n = n ? n : m
matrix e = vec(mshape(seq(1, m*n), m, n)')
return post ? A[,e] : A[e,]
```

```
# Generates truncated normal random values.
# Set 'below' and/or 'above' to NA to skip.
scalar l = ok(below) ? cnorm((below - m)/sigma) : 0
scalar r = ok(above) ? cnorm((above - m)/sigma) : 1
matrix u = l + (r-l) .* muniform(n,1)
return invcdf(z, u) .* sigma + m
```

```
series DOK = diff(ok(y))
series Beg = DOK == -1
series End = DOK(+1) == 1
series splen = 0
splen = DOK == -1 ? 1 : DOK==1 ? 0 : (splen(-1) == 0 ? 0 : splen(-1) + 1)
series y0 = NA
series y0 = Beg ? y(-1) : y0(-1)
series y1 = NA
series y1 = End ? y(+1) : y1(-1)
set skip_missing off
matrix A = {y, y0, y1, splen}
scalar t = lastobs(y)
loop while t>firstobs(y) --quiet
if ok(End[t]) && (End[t] == 1)
scalar l = A[t, 4]
dy = (A[t,3] - A[t,2]) / (l + 1)
patch = A[t,2] + dy * seq(1,l)'
A[t-l+1:t,1] = patch
t -= l
else
t--
endif
endloop
return A[,1]
```