text:
export printuln, usprint, save2fig, save3fig

@doc raw"""
# Function call

`printuln(r...; label="")`

# Description

This function prints a real or complex variable in an optionally specified unit

# Variables

`r[1]` Text string to printed, indicating the variable name or description; a
maximum of 16 characters is allowed, in case the optional variable `label` is
NOT specified; in  case of specifying the variable `label`, the maximum length
of `r[1]` may not exceed 12 characters

`r[2]` Variable to be printed; the variable may be real or complex, units
assigned from the module Unitful are allowed; vectors of variables are also
allowed

`r[3]` Optional target unit of `r[2]` to be displayed

`label` Optional four character label structuring the output of variables, e.g.,
`label = "(a)"`

# Examples

```julia
julia> using Unitful, Unitful.DefaultSymbols, ElectricalEngineering
julia> U1 = 300V+j*400V
julia> printuln("U1", U1, kV)
              U1 = 0.3 kV + j 0.4 kV
                 = 0.5 kV ∠ 53.1301°
julia> printuln("real(U1)", real(U1), kV)
        real(U1) = 0.3 kV
julia> printuln("U1", U1, V, label="(a)")
(a)           U1 = 300 V + j 400 V
                 = 500 V ∠ 53.1301°
```
"""
function printuln(r...; label="")
    # r[1] = 1st argument = string to be printed
    # r[2] = 2nd argument = quantity to shown
    # r[3] = 3rd argument = optional output unit
    # If unit is specified as preferred 3rd argument, then convert number
    # to this argument
    for k in collect(1:size(r[2], 1))
        if length(r) >= 3
            # Convert numerical value to specified units
            # Change printed unit
            num = uconvert(r[3], r[2][k])
        else
            num = r[2][k]
        end
        numunit = unit(num)

        # Real and imaginary part
        numre = ustrip(real(num))
        numim = ustrip(imag(num))
        # Magnitude and angle
        numabs = abs(ustrip(num))
        numargdeg = angle(ustrip(num))*180/pi
        # Determine unit associated to numeric value and consider
        # special spelling of units

        # Print text string (1st argument)
        if size(r[2], 1) == 1
          # Second argument (string) is a scalar value
          if label == ""
              @printf("%16s = ", r[1])
          else
              @printf("%-4s%12s = ",label, r[1])
          end
        else
          # Second argument (string) is a vector
          if label == ""
              @printf("%16s = ", r[1]*"["*@sprintf("%i", k)*"]")
          else
              if k == 1
                  # Print label only in first entry
                  @printf("%-4s%12s = ", label, r[1]*"["*@sprintf("%i", k)*"]")
              else
                  # Omit label on sencond and following entries
                  @printf("%16s = ", r[1]*"["*@sprintf("%i", k)*"]")
              end
          end
        end
        @printf("%g %s", numre, numunit)
        if numim > 0
            @printf(" + j*%g %s\n", abs(numim), numunit)
            @printf("%16s = %g %s ∠ %g°\n","",numabs, numunit, numargdeg)
        elseif ustrip(imag(num)) < 0
            @printf(" - j*%g %s\n",abs(numim), numunit)
            @printf("%16s = %g %s ∠ %g°\n","",numabs, numunit, numargdeg)
        else
            @printf("\n")
        end
    end
end

"""
# Function call

`usprint(u,U)`

# Description

Returns a string containing the number with six digits plus unit without space
in between, such that the string output can be copied and pasted.

# Variables

`u` Variable to be processed

`U` Optional target unit; if not provided, the actual unit is used

# Examples
```julia
julia> I1 = 1mA
julia> usprint(I1)
1mA
julia> usprint(I1, A)
0.001A
"""
function usprint(u, U=unit(u))
    q = uconvert(U, u)
    return @sprintf("%g%s", Float64(ustrip(q)), U)
end

"""
# Function call

`save2fig(fileName, subDir="."; dpi=300, tight=true, crop=false)`

# Description

Save the actual figure in the following three file formats
- `png` Portable network graphics in subdirectory png/
- `pdf` Encapsulated postscript in subdirectory pdf/
These two graphics format are relevant when processing figures in LaTeX,
LyX, LibreOffice or other word processors

# Variables

`fileName` String of the filename, extended by `.png`, and `.pdf`; the
file names are stored in sub directories `png/`, and `pdf/`

`subDir` String of sub directory in which the adjacent sub directories  `png/`,
and `pdf/` are located; default = `.` (local work directory)

`dpi` Resolution of the stored PNG file; default = 300 (dpi)

`tight` Additional margin around the figure is removed, except for a 2% margin;
default = `true`

`crop` Crop all 'white' space around the figure; this feature requires the
installation of the following software
- imagemagick to be applied to `png` files,
  see https://www.imagemagick.org/script/index.php on Windows
  or apply `sudo apt install imagemagick` on Linux (Mint or Ubuntu)
- pdfcrop to be applied to `pdf` files, see
  https://www.ctan.org/pkg/pdfcrop on Windows
  or apply `sudo apt install texlive-extra-utils` on Linux (Mint or Ubuntu)

# Examples

```julia
julia> save2fig("phasordiagram")
julia> save2fig("phasordiagram_crop", crop=true)
```
"""
function save2fig(fileName, subDir="."; dpi=300, tight=true, crop=false)
    # Store PNG file
    mkpath(subDir*"/png")
    if tight == false
        savefig(subDir*"/png/"*fileName*".png", dpi=dpi)
    else
        savefig(subDir*"/png/"*fileName*".png", dpi=dpi,
            bbox_inches = "tight", pad_inches=0.02)
    end
    if crop==true
        try
            arg = `$subDir/png/$fileName.png`
            status = run(`convert $arg -trim $arg`);
        catch err
            error("module ElectricalEngineering: function save3fig: Binary file not found: convert
    The software convert (imagemagick) may not be installed or the path may not
    be specified or the software may have caused a runtime error
    To install imagemagick, see https://www.imagemagick.org/script/index.php on
    Windows or apply `sudo apt install imagemagick` on Linux (Mint or Ubuntu)")
        end
    end

    # Store PDF file
    mkpath(subDir*"/pdf")
    if tight == false
        savefig(subDir*"/pdf/"*fileName*".pdf")
    else
        savefig(subDir*"/pdf/"*fileName*".pdf",
            bbox_inches="tight", pad_inches=0.02)
    end
    if crop == true
        try
            arg = `$subDir/pdf/$fileName.pdf`
            status=run(`pdfcrop $arg $arg`);
        catch err
            error("module ElectricalEngineering: function save3fig: Binary file not found: dpfcrop
    The software pdfcrop may not be installed or the path may not be specified
    or the software may have caused a runtime error
    To install pdfcrop, see https://www.ctan.org/pkg/pdfcrop on Windows
    or apply `sudo apt install texlive-extra-utils` on Linux (Mint or Ubuntu)")
        end
    end
end

"""
# Function call

`save3fig(fileName, subDir="."; dpi=300, tight=true, crop=false)`

# Description

Save the actual figure in the following three file formats
- `png` Portable network graphics in subdirectory png/
- `eps` Encapsulated postscript in subdirectory eps/
- `pdf` Encapsulated postscript in subdirectory pdf/
These three graphics format are relevant when processing figures in LaTeX,
LyX, LibreOffice or other word processors

# Variables

`fileName` String of the filename, extended by `.png`, `.eps` and `.pdf`; the
file names are stored in sub directories `png/`, `eps/` and `pdf/`

`subDir` String of sub directory in which the adjacent sub directories  `png/`,
`eps/` and `pdf/` are located; default = `.` (local work directory)

`dpi` Resolution of the stored PNG file; default = 300 (dpi)

`tight` Additional margin around the figure is removed, except for a 2% margin;
default = `true`

`crop` Crop all 'white' space around the figure; this feature requires the
installation of the following software
- imagemagick to be applied to `png` files,
  see https://www.imagemagick.org/script/index.php on Windows
  or apply `sudo apt install imagemagick` on Linux (Mint or Ubuntu)
- epstool to be applied to `eps` files, see
  http://pages.cs.wisc.edu/~ghost/gsview/epstool.htm on Windows
  or apply `sudo apt install epstool` on Linux (Mint or Ubuntu)
- pdfcrop to be applied to `pdf` files, see
  https://www.ctan.org/pkg/pdfcrop on Windows
  or apply `sudo apt install texlive-extra-utils` on Linux (Mint or Ubuntu)

# Examples

```julia
julia> save3fig("phasordiagram")
julia> save3fig("phasordiagram_crop", crop=true)
```
"""
function save3fig(fileName, subDir="."; dpi=300, tight=true, crop=false)
    save2fig(fileName, subDir = subDir, dpi = dpi, tight = tight, crop = crop)

    # Store EPS file
    mkpath(subDir*"/eps")
    if tight == false
        savefig(subDir*"/eps/"*fileName*".eps")
    else
        savefig(subDir*"/eps/"*fileName*".eps",
            bbox_inches = "tight", pad_inches=0.02)
    end
    if crop==true
        try
            arg = `$subDir/eps/$fileName.eps`
            argtemp = `$subDir/eps/$fileName.eps.temp`
            status = run(`epstool --copy --bbox $arg $argtemp`);
            mv(subDir*"/eps/"*fileName*".eps.temp",
                subDir*"/eps/"*fileName*".eps",
                remove_destination=true)
        catch err
            error("module ElectricalEngineering: function save3fig: Binary file not found: epstool
    The software epstool may not be installed or the path may not be specified
    or the software may have caused a runtime error
    To install epstool, see http://pages.cs.wisc.edu/~ghost/gsview/epstool.htm
    on Windows or apply `sudo apt install epstool` on Linux (Mint or Ubuntu)")
        end
    end
end


text:
@testset "can_avx" begin

  @test LoopVectorization.ArrayInterface.can_avx(log)
  @test LoopVectorization.ArrayInterface.can_avx(log1p)
  @test LoopVectorization.ArrayInterface.can_avx(exp)
  @test LoopVectorization.ArrayInterface.can_avx(+)
  @test LoopVectorization.ArrayInterface.can_avx(-)
  @test LoopVectorization.ArrayInterface.can_avx(Base.FastMath.add_fast)
  @test LoopVectorization.ArrayInterface.can_avx(/)
  @test LoopVectorization.ArrayInterface.can_avx(sqrt)
  @test LoopVectorization.ArrayInterface.can_avx(tanh_fast)
  @test LoopVectorization.ArrayInterface.can_avx(sigmoid_fast)
  @test LoopVectorization.ArrayInterface.can_avx(LoopVectorization.relu)
  @test !LoopVectorization.ArrayInterface.can_avx(clenshaw)
  @test !LoopVectorization.ArrayInterface.can_avx(println)

end

ext: jl


text:
## QUESTION 2

include("md.jl");

# Pick interaction distance
δ = 0.1

# Create a new simulation object
qtherm = MolecularDynamicsTrial(
    1000,               # number of bodies
    0.75,               # particle density
    1.069,              # desired temperature; also use 1.304
    0.032,              # time step dx, in units of τ
    2000,              # number of iterations
    3.405e-10,          # length units, from Verlet, in m
    119.8,              # energy scale for potential, in Kelvin
    39.948 / 6.0221e26  # argon atomic mass, in kg
);

# Place the particles
mdplace!(qtherm);

# Run for 2000 steps
metropolis!(qtherm, 1999, δ)

# Save the trial
print("Done thermalizing. Saving new starting position...")
save("../../swp/1069_met_therm.mdv", qtherm)

# Take the last step as the start step
qstart = subtrial(qtherm, 2000, 2000)
save("../../swp/1069_met_start.mdv", qstart)
println("done.")


text:
using ConstrainedDynamics

Mechanism("urdf_test.urdf", floating=false)
@test true

Mechanism("urdf_test.urdf", floating=true)
@test true


text:
const severity_hierarchy = ["debug", "info", "warn", "error"]

const message_severities = [
    DEBUG_UTILS_MESSAGE_SEVERITY_VERBOSE_BIT_EXT,
    DEBUG_UTILS_MESSAGE_SEVERITY_INFO_BIT_EXT,
    DEBUG_UTILS_MESSAGE_SEVERITY_WARNING_BIT_EXT,
    DEBUG_UTILS_MESSAGE_SEVERITY_ERROR_BIT_EXT,
]

const message_types = Dict(
    "general" => DEBUG_UTILS_MESSAGE_TYPE_GENERAL_BIT_EXT,
    "validation" => DEBUG_UTILS_MESSAGE_TYPE_VALIDATION_BIT_EXT,
    "performance" => DEBUG_UTILS_MESSAGE_TYPE_PERFORMANCE_BIT_EXT,
)

const message_types_r = reverse_dict(message_types)

function vk_log(message_severity, str)
    log = @match message_severity begin
        &DEBUG_UTILS_MESSAGE_SEVERITY_VERBOSE_BIT_EXT => "\e[34;1;1m[ Debug"
        &DEBUG_UTILS_MESSAGE_SEVERITY_INFO_BIT_EXT => "\e[36;1;1m[ Info"
        &DEBUG_UTILS_MESSAGE_SEVERITY_WARNING_BIT_EXT => "\e[33;1;1m[ Warning"
        &DEBUG_UTILS_MESSAGE_SEVERITY_ERROR_BIT_EXT => "\e[31;1;1m[ Error"
        _ => error("Unknown message severity $message_severity")
    end
    ccall(:jl_safe_printf, Cvoid, (Cstring,), string(log, ": \e[m", str, '\n'))
end

"""
Default callback for debugging with [`DebugUtilsMessengerEXT`](@ref).
"""
function default_debug_callback(message_severity, message_type, callback_data_ptr, user_data_ptr)
    callback_data_ptr == C_NULL && return UInt32(0)
    callback_data = unsafe_load(callback_data_ptr)
    message = unsafe_string(callback_data.pMessage)
    message_severity = DebugUtilsMessageSeverityFlagEXT(UInt32(message_severity))
    message_type = DebugUtilsMessageTypeFlagEXT(UInt32(message_type))
    if !startswith(message, "Device Extension: VK")
        id_name = unsafe_string(callback_data.pMessageIdName)
        vk_log(message_severity, "$(uppercase(message_types_r[message_type])) ($id_name): $message")
    end
    UInt32(0)
end

"""
Register a user-defined callback and return the corresponding messenger.

A default function [`default_debug_callback`](@ref) can be converted to a function pointer to use as a callback.

!!! warning
    `callback` must be a function pointer of type `Ptr{Nothing}` obtained from a `callback_f` function as follows:
    `callback = @cfunction(callback_f, UInt32, (VkDebugUtilsMessageSeverityFlagBitsEXT, VkDebugUtilsMessageTypeFlagBitsEXT, Ptr{vk.VkDebugUtilsMessengerCallbackDataEXT}, Ptr{Cvoid}))`
    with `callback_f` a Julia function with a signature matching the `@cfunction` call.
"""
function DebugUtilsMessengerEXT(instance::Instance, callback::Ptr{Nothing}; severity = "info", types = ["general", "validation", "performance"])
    index = findfirst(severity_hierarchy .== severity)
    severity_bits = message_severities[index:end]
    type_bits = [message_types[key] for key in types]
    DebugUtilsMessengerEXT(instance, |(severity_bits...), |(type_bits...), callback, function_pointer(instance, "vkCreateDebugUtilsMessengerEXT"), function_pointer(instance, "vkDestroyDebugUtilsMessengerEXT"))
end


text:
"""
    Benchmarking for PDIPS.jl.

Note that benchmarking times do not include problem reformulation.
"""

using BenchmarkTools
using MatrixDepot
const MD = MatrixDepot
using Printf

using PDIPS

const NETLIB_PROBLEMS = [
    "LPnetlib/lp_25fv47",
    "LPnetlib/lp_adlittle",
    "LPnetlib/lp_afiro",
    "LPnetlib/lp_agg",
    "LPnetlib/lp_brandy",
    "LPnetlib/lpi_chemcom",
    "LPnetlib/lp_fit1d",
    "LPnetlib/lp_ganges",
    "LPnetlib/lp_stocfor1"
]

function load_descriptors!(mds::Vector{MD.MatrixDescriptor},
                        names::Vector{String};
                        verbose::Bool = true)
    np = length(names)

    for i in Base.OneTo(np)
        verbose && @printf("Loading %s...\n", names[i])
        mds[i] = mdopen(names[i])
    end

    nothing
end

function reformat(vector::Vector{Float64})
    n = length(vector)
    _vec = Vector{Float64}(undef, n)
    for i in Base.OneTo(n)
        _vec[i] = vector[i] == 1.0e308 ? Inf : vector[i]
    end

    _vec
end

# Problems from Netlib
PROBLEMS = Vector{MD.MatrixDescriptor}(undef, length(NETLIB_PROBLEMS))
load_descriptors!(PROBLEMS, NETLIB_PROBLEMS)

#==============
    25fv47
==============#
p = PROBLEMS[1]

# load problem
lp = Problem{Float64}()
load_problem!(lp, p.A, vec(p.b), vec(p.c), vec(p.lo), reformat(vec(p.hi)))
bm = @benchmarkable solve(lp)
tune!(bm)
run(bm)

solve(lp)

#===============
    adlittle
===============#
p = PROBLEMS[2]

# load problem
lp = Problem{Float64}()
load_problem!(lp, p.A, vec(p.b), vec(p.c), vec(p.lo), reformat(vec(p.hi)))
bm = @benchmarkable solve(lp)
tune!(bm)
run(bm)

#============
    afiro
============#
p = PROBLEMS[3]

# load problem
lp = Problem{Float64}()
load_problem!(lp, p.A, vec(p.b), vec(p.c), vec(p.lo), reformat(vec(p.hi)))
bm = @benchmarkable solve(lp)
tune!(bm)
run(bm)

#==========
    agg
==========#
p = PROBLEMS[4]

# load problem
lp = Problem{Float64}()
load_problem!(lp, p.A, vec(p.b), vec(p.c), vec(p.lo), reformat(vec(p.hi)))
bm = @benchmarkable solve(lp)
tune!(bm)
run(bm)

#=============
    brandy
=============#
p = PROBLEMS[5]

# load problem
lp = Problem{Float64}()
load_problem!(lp, p.A, vec(p.b), vec(p.c), vec(p.lo), reformat(vec(p.hi)))
bm = @benchmarkable solve(lp)
tune!(bm)
run(bm)

#==============
    chemcom
==============#
p = PROBLEMS[6]

# load problem
lp = Problem{Float64}()
load_problem!(lp, p.A, vec(p.b), vec(p.c), vec(p.lo), reformat(vec(p.hi)))
bm = @benchmarkable solve(lp)
tune!(bm)
run(bm)

#=============
    fit1d
=============#
p = PROBLEMS[7]

# load problem
lp = Problem{Float64}()
load_problem!(lp, p.A, vec(p.b), vec(p.c), vec(p.lo), reformat(vec(p.hi)))
bm = @benchmarkable solve(lp)
tune!(bm)
run(bm)

#=============
    ganges
=============#
p = PROBLEMS[8]

# load problem
lp = Problem{Float64}()
load_problem!(lp, p.A, vec(p.b), vec(p.c), vec(p.lo), reformat(vec(p.hi)))
bm = @benchmarkable solve(lp)
tune!(bm)
run(bm)

#===============
    stocfor1
===============#
p = PROBLEMS[9]

# load problem
lp = Problem{Float64}()
load_problem!(lp, p.A, vec(p.b), vec(p.c), vec(p.lo), reformat(vec(p.hi)))
bm = @benchmarkable solve(lp)
tune!(bm)
run(bm)


text:
# module Lupin

# code from julia/base/namedtuple.jl
#           julia/base/abstractdict.jl

if VERSION >= v"1.5.0-DEV.182"
import Base: mergewith
end
using Base: Callable, merge_names, merge_types, sym_in

function mergewith(combine::Callable, nt::NamedTuple, others::NamedTuple...)
    foldl(mergewith(combine), others; init = nt)
end

function mergewith(combine::Callable, a::NamedTuple{an}, b::NamedTuple{bn}) where {an, bn}
    names = merge_names(an, bn)
    types = collect(fieldtypes(merge_types(names, typeof(a), typeof(b))))
    values = @inbounds map(enumerate(names)) do (idx, n)
            if sym_in(n, bn)
                if sym_in(n, an)
                    value = combine(getfield(a, n), getfield(b, n))
                    typ = types[idx]
                    if typ !== typeof(value)
                        types[idx] = promote_type(typ, typeof(value))
                    end
                    value
                else
                    getfield(b, n)
                end
            else
                getfield(a, n)
            end
        end
    NamedTuple{names, Tuple{types...}}(values)
end

if VERSION < v"1.5.0-DEV.182"
    mergewith(combine::Callable) = (args::NamedTuple...) -> mergewith(combine, args...)
    export mergewith
end

# module Lupin


text:
module Opcode

export AABBTree, intersect, intersect!, deform!, raycast, raycast!

using CxxWrap

const depsfile = joinpath(dirname(dirname(@__FILE__)), "deps", "deps.jl")
if !isfile(depsfile)
    error("$depsfile not found, Opcode did not build properly")
end
include(depsfile)
const jlcxx_path = _l_jlcxx

# Wrap the functions defined in C++
wrap_module(Opcode._l_opcode, Opcode)

const AABBTree = Opcode.ModelRef

function AABBTree(verts::Array{<:Real,2}, faces::Array{<:Integer,2})
    size(verts,1)==3 && size(faces,1)==3 || error("require verts : 3 x n and faces : 3 x m")
    AABBTree(Opcode.aabb_tree_create(verts,size(verts,2),faces,size(faces,2)))
end

# cast first n rays in (origins,dirs)
function raycast(tree::AABBTree, origins::Array{<:Real,2}, dirs::Array{<:Real,2}, n::Integer = size(origins,2);
                 mode::Int=0)
    size(origins,1)==size(dirs,1)==3 &&
    n <= min(size(origins,2), size(dirs,2)) || error("require origins : 3 x n and dirs : 3 x n")
    if mode==0
        hits = Vector{Int}()
    	dists = Vector{Float64}();
    	faceinds = Vector{Int}();
    	barys = Vector{Float64}();
    	Opcode.aabb_tree_intersect(tree, origins, dirs, hits, dists, faceinds, barys, n)
    	(hits,dists,faceinds,reshape(barys,2,length(hits)))
    else
        hits = Vector{Int}(n)
        dists = Vector{Float64}(n);
        faceinds = Vector{Int}(n);
        barys = Matrix{Float64}(2,n);
        Opcode.aabb_tree_intersect_mut(tree, origins, dirs, hits, dists, faceinds, barys, n)
        (hits,dists,faceinds,reshape(barys,2,length(hits)))
    end
end

function raycast!(tree::AABBTree, origins::Array{<:Real,2}, dirs::Array{<:Real,2},
                    hits::Vector{<:Integer}, dists::Vector{<:Real}, faceinds::Vector{<:Integer},
		    barys::Array{<:Real,2}, n :: Integer = size(origins,2))
    size(origins,1)==size(dirs,1)==3 && size(barys,1)==2 &&
    n <= min(size(origins,2),size(dirs,2),length(hits),length(dists),length(faceinds),size(barys,2)) ||
    error("input dimensions")
    Opcode.aabb_tree_intersect_mut(tree, origins, dirs, hits, dists, faceinds, barys, n)
    tree
end

function deform!(tree::AABBTree, verts::Array{<:Real,2})
    size(verts,1)==3 || error("require verts : 3 x n")
    Opcode.aabb_tree_deform(tree, verts, size(verts,2))
    tree
end

function bary2eucl!(x::Array{<:Real,2}, bary::Array{<:Real,2}, faceinds::Vector{<:Integer},
                   verts::Array{<:Real,2}, faces::Array{<:Integer})
    n = size(bary,2)
    size(bary,1)==2 || error("bary : 2 x n")
    @simd for i = 1:n
    	b2 = bary[1,i]
	b3 = bary[2,i]
        b1 = 1 - b2 - b3
        v1,v2,v3 = (faces[1,faceinds[i]], faces[2,faceinds[i]], faces[3,faceinds[i]])
        v1x,v1y,v1z = (verts[1,v1], verts[2,v1], verts[3,v1])
        v2x,v2y,v2z = (verts[1,v2], verts[2,v2], verts[3,v2])
        v3x,v3y,v3z = (verts[1,v3], verts[2,v3], verts[3,v3])
        x[1,i] = b1*v1x + b2*v2x + b3*v3x
        x[2,i] = b1*v1y + b2*v2y + b3*v3y
        x[3,i] = b1*v1z + b2*v2z + b3*v3z
    end
    x
end
function bary2eucl!(bary::Array{T,2}, faceinds::Vector{<:Integer},
                    verts::Array{<:Real,2}, faces::Array{<:Integer}) where {T<:Real}
    x = Array{T,2}(3,n)
    bary2eucl(x,bary,faceinds,verts,faces)
end

end # module


text:
export zeropad
using AutoGrad

# cut-and-paste method: see Knet/prof/zeropad for alternatives and timing

function zeropad(x, p; y=nothing, fillzero=true)
    p = padtuple(p, size(x))
    n = ndims(x)
    d = ntuple(i->(i < n-1 ? size(x,i)+sum(p[i]) : size(x,i)), n)
    c = ntuple(i->(i < n-1 ? (p[i][1]+1:p[i][1]+size(x,i)) : Colon()), n)
    y === nothing ? y = similar(x,d) : @assert typeof(y)===typeof(x) && size(y)==d
    fillzero && fill!(y, 0)
    y[c...] .= x
    return y
end

# Gradients

function ∇zeropad(dy, p)
    p = padtuple(p, size(dy))
    n = ndims(dy)
    c = ntuple(i->(i < n-1 ? (p[i][1]+1:size(dy,i)-p[i][2]) : Colon()), n)
    dy[c...]
end

@primitive1 zeropad(x,p...;o...),dy,y  ∇zeropad(dy,p...)
@primitive1 ∇zeropad(dy,p...),ddx,dx  zeropad(ddx,p...)

# padtuple turns Int/Vector/Tuple to ((a,b),(c,d)) format

# Check if the right ((a,b),(c,d)) format
padtuple(p::NTuple{P,NTuple{2,Int}}, d::Dims{D}) where {P,D} = (@assert P == D-2; p)

# Main conversions
padtuple(p::Int, d::Dims{D}) where {D} = ntuple(i->(p,p), D-2)
padtuple(p::Vector{Int}, d::Dims{D}) where {D} = (@assert length(p) == D-2; ntuple(i->(p[i],p[i]), D-2))
padtuple(p::NTuple{P,Int}, d::Dims{D}) where {P,D} = (@assert P == D-2; ntuple(i->(p[i],p[i]), P))

# Handle other integer types
padtuple(p::Integer, d) = padtuple(Int(p), d)
padtuple(p::Vector{<:Integer}, d) = padtuple(Int.(p), d)
padtuple(p::NTuple{P,<:Integer}, d) where {P} = padtuple(Int.(p), d)
padtuple(p::NTuple{P,NTuple{2,<:Integer}}, d) where {P} = padtuple(ntuple(i->Int.(p[i]),P), d)

# symmetric padding can be handled by the convolution
symmetric(p::Integer) = true
symmetric(p::Vector{<:Integer}) = true
symmetric(p::NTuple{P,<:Integer}) where {P} = true
symmetric(p::NTuple{P,NTuple{2,<:Integer}}) where {P} = all(x[1]==x[2] for x in p)


text:
using Distributions

mutable struct Simulation{N}
    heatmap::Array{Float64, N}
    stats::NamedTuple{(:mean, :std)}
end

function simulate(distribution, dims, n)
    tp = ntuple(i -> n, dims)
    heatmap = rand(distribution, tp...)
    return Simulation{dims}(heatmap, (mean = mean(heatmap), std = std(heatmap)))
end

sim = simulate(Normal(), 2, 1000);

#=
julia> sim = simulate(Normal(), 2, 1000);

julia> sim.heatmap
1000×1000 Array{Float64,2}:
  1.22768    -0.316595    1.53702     …   0.944001    1.14345
 -1.09027    -0.553394    1.23286     …   0.779135   -0.966926
  1.52175     0.297313    1.02013         0.6078     -0.299004
  ⋮                                   ⋱
  0.912278   -0.219911    0.513086       -1.37209     0.126663
  1.05601     0.650547   -0.490138       -0.321995    0.455638
 -0.764425    0.175759   -0.245518        0.724038   -0.637692

julia> sim.stats
(mean = 0.0003489176881672829, std = 1.0002976677059898)
=#

# accessors - getters
# they are just forwarding to the properties

get_heatmap(s::Simulation) = s.heatmap
get_stats(s::Simulation) = s.stats

# alternative naming convention
# note the extension of Base.size generic function
heatmap(s::Simulation) = s.heatmap
stats(s::Simulation) = s.stats

# accessors - setters
# we only want a setter for `points` but not `size`
# we can also add validation.

function heatmap!(
        s::Simulation{N},
        new_heatmap::AbstractArray{Float64, N}) where {N}
    s.heatmap = new_heatmap
    s.stats = (mean = mean(new_heatmap), std = std(new_heatmap))
    return nothing
end

# add validation
function heatmap!(
            s::Simulation{N},
            new_heatmap::AbstractArray{Float64, N}) where {N}
    if length(unique(size(new_heatmap))) != 1
        error("All dimensions must have same size")
    end
    s.heatmap = new_heatmap
    s.stats = (mean = mean(new_heatmap), std = std(new_heatmap))
    return nothing
end

#=
julia> heatmap!(sim, rand(10, 3))
ERROR: All dimensions must have same size
=#

# ------------------------------------------------------------------------
# marking fields private

mutable struct Simulation{N}
    _heatmap::Array{Float64, N}
    _stats::NamedTuple{(:mean, :std)}
end


text:
module DecisionMaking

# Write your package code here.

end


text:
##############################################################################
#
# test_2d_vortex_particles.jl
#
# Check 2D vortex particles behave as expected.
#
# Copyright 2020 HJA Bird
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
# IN THE SOFTWARE.
##############################################################################

function two_particles_basic_case_2d()
	@testset "2D 2 particle interaction" begin
		pl0 = [0,0]
		ply = [0,1]
		plx = [1,0]

		# Singular kernel ---------------------------------------------------------------------------------------
		kernel = singular_regularisation()
		@test isapprox([1,0]./(2*pi), particle_induced_velocity(pl0, 1, ply, kernel, 1.); rtol=1e-5)
		@test isapprox([0,-1]./(2*pi), particle_induced_velocity(pl0, 1, plx, kernel, 1.); rtol=1e-5)
		@test isapprox([2,0]./(2*pi), particle_induced_velocity(pl0, 2, ply, kernel, 1.); rtol=1e-5)
        @test isapprox([0.5,0]./(2*pi), particle_induced_velocity(pl0, 1, 2*ply, kernel, 1.); rtol=1e-5)
	end
end

function test_correct_output_matrix_dims_2d_vortex_particles()
    @testset "2D particle correct output dimensions" begin
        n1 = 31
        n2 = 47
        kernel = singular_regularisation()
        @test (2,) == size(particle_induced_velocity(rand(2), rand(), rand(2), kernel, 1.))
        @test (2,) == size(particle_induced_velocity(rand(n1, 2), rand(n1), rand(2), kernel, 1.))
        @test (n1, 2) == size(particle_induced_velocity(rand(n1, 2), rand(n1), rand(n1, 2), kernel, 1.))
        @test (n2, 2) == size(particle_induced_velocity(rand(n1, 2), rand(n1), rand(n2, 2), kernel, 1.))
        @test (n1, 2) == size(particle_induced_velocity(rand(2), rand(), rand(n1, 2), kernel, 1.))
        kernel = gaussian_regularisation()
        @test Float32 == typeof(particle_visc_induced_dvort(
            rand(2), rand(), rand(), rand(2), rand(), rand(), kernel, 1., 1))
        @test Float32 == typeof(particle_visc_induced_dvort(
            rand(n1, 2), rand(n1), rand(n1), rand(2), rand(), rand(), kernel, 1., 1))
        @test (n1,) == size(particle_visc_induced_dvort(
            rand(n1, 2), rand(n1), rand(n1), rand(n1, 2), rand(n1), rand(n1), kernel, 1., 1))
        @test (n2,) == size(particle_visc_induced_dvort(
            rand(n1, 2), rand(n1), rand(n1), rand(n2, 2), rand(n2), rand(n2), kernel, 1., 1))
        @test (n1,) == size(particle_visc_induced_dvort(
            rand(2), rand(), rand(), rand(n1, 2), rand(n1), rand(n1), kernel, 1., 1))
    end
end

test_correct_output_matrix_dims_2d_vortex_particles()
two_particles_basic_case_2d()


text:
using PlotThemes
using Test, PlotUtils

@test in(:sand, keys(PlotThemes._themes))


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""Spec to control the desired behavior of daemon set rolling update.

    IoK8sApiAppsV1RollingUpdateDaemonSet(;
        maxUnavailable=nothing,
    )

    - maxUnavailable::IoK8sApimachineryPkgUtilIntstrIntOrString : The maximum number of DaemonSet pods that can be unavailable during the update. Value can be an absolute number (ex: 5) or a percentage of total number of DaemonSet pods at the start of the update (ex: 10%). Absolute number is calculated from percentage by rounding up. This cannot be 0. Default value is 1. Example: when this is set to 30%, at most 30% of the total number of nodes that should be running the daemon pod (i.e. status.desiredNumberScheduled) can have their pods stopped for an update at any given time. The update starts by stopping at most 30% of those DaemonSet pods and then brings up new DaemonSet pods in their place. Once the new pods are available, it then proceeds onto other DaemonSet pods, thus ensuring that at least 70% of original number of DaemonSet pods are available at all times during the update.
"""
mutable struct IoK8sApiAppsV1RollingUpdateDaemonSet <: SwaggerModel
    maxUnavailable::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgUtilIntstrIntOrString } # spec name: maxUnavailable

    function IoK8sApiAppsV1RollingUpdateDaemonSet(;maxUnavailable=nothing)
        o = new()
        validate_property(IoK8sApiAppsV1RollingUpdateDaemonSet, Symbol("maxUnavailable"), maxUnavailable)
        setfield!(o, Symbol("maxUnavailable"), maxUnavailable)
        o
    end
end # type IoK8sApiAppsV1RollingUpdateDaemonSet

const _property_map_IoK8sApiAppsV1RollingUpdateDaemonSet = Dict{Symbol,Symbol}(Symbol("maxUnavailable")=>Symbol("maxUnavailable"))
const _property_types_IoK8sApiAppsV1RollingUpdateDaemonSet = Dict{Symbol,String}(Symbol("maxUnavailable")=>"IoK8sApimachineryPkgUtilIntstrIntOrString")
Base.propertynames(::Type{ IoK8sApiAppsV1RollingUpdateDaemonSet }) = collect(keys(_property_map_IoK8sApiAppsV1RollingUpdateDaemonSet))
Swagger.property_type(::Type{ IoK8sApiAppsV1RollingUpdateDaemonSet }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAppsV1RollingUpdateDaemonSet[name]))}
Swagger.field_name(::Type{ IoK8sApiAppsV1RollingUpdateDaemonSet }, property_name::Symbol) =  _property_map_IoK8sApiAppsV1RollingUpdateDaemonSet[property_name]

function check_required(o::IoK8sApiAppsV1RollingUpdateDaemonSet)
    true
end

function validate_property(::Type{ IoK8sApiAppsV1RollingUpdateDaemonSet }, name::Symbol, val)
end


text:
# Write function to automate plotting combined diagnostics plots
# - Plots all three standard plots for assessing linear model fits

function plot_model_check(;model, df, y)

    ###################################################
    # - Setup and defaults
    ###################################################

    # Set Gadfly theme
    default_theme = Gadfly.Theme(
        panel_fill = nothing,
        highlight_width = 0Gadfly.pt,
        point_size = 0.8Gadfly.mm,
        key_position = :inside,
        #grid_line_width = [0 Measures.mm],
        panel_stroke = Colors.colorant"black",
        background_color = Colors.colorant"white")
    Gadfly.push_theme(default_theme)

    # Set default plot size
    Gadfly.set_default_plot_size(24Gadfly.cm, 12Gadfly.cm)

    # Make individual diagnostics plots
    # - These plots are the equivalent of calling plot(...) on a
    #   linear model (e.g. lm, ANOVA) in R.
    p1 = ModelCheck.plot_qqplot(model = model, df = df, y = y)
    p2 = ModelCheck.plot_residual_fitted(model = model, df = df, y = y)
    p3 = ModelCheck.plot_scale_location(model = model, df = df, y = y)

    # Plot all three plots in a single graph
    all_plots = Gadfly.hstack(p1, p2, p3)

    # Return plots
    return(all_plots)

end


text:
module Runoff

using LinearAlgebra
using DocStringExtensions

using ...VariableTemplates
using ...Land: SoilModel, pressure_head, hydraulic_conductivity, get_temperature

export AbstractPrecipModel,
    DrivenConstantPrecip,
    AbstractSurfaceRunoffModel,
    NoRunoff,
    compute_surface_grad_bc,
    CoarseGridRunoff

"""
    AbstractPrecipModel{FT <: AbstractFloat}
"""
abstract type AbstractPrecipModel{FT <: AbstractFloat} end

"""
    DrivenConstantPrecip{FT, F} <: AbstractPrecipModel{FT}

Instance of a precipitation distribution where the precipication value
is constant across the domain. However, this value can change in time.

# Fields
$(DocStringExtensions.FIELDS)
"""
struct DrivenConstantPrecip{FT, F} <: AbstractPrecipModel{FT}
    "Mean precipitation in grid"
    mp::F
    function DrivenConstantPrecip{FT}(mp::F) where {FT, F}
        new{FT, F}(mp)
    end
end

function (dcp::DrivenConstantPrecip{FT})(t::Real) where {FT}
    return FT(dcp.mp(t))
end

"""
    AbstractSurfaceRunoffModel

Abstract type for different surface runoff models. Currently, only
`NoRunoff` is supported.
"""
abstract type AbstractSurfaceRunoffModel end

"""
    NoRunoff <: AbstractSurfaceRunoffModel

Chosen when no runoff is to be modeled.
"""
struct NoRunoff <: AbstractSurfaceRunoffModel end

"""
    CoarseGridRunoff{FT} <: AbstractSurfaceRunoffModel

Chosen when no subgrid effects are to be modeled.
"""
struct CoarseGridRunoff{FT} <: AbstractSurfaceRunoffModel
    "Mean vertical resolution at the surface"
    Δz::FT
end

"""
    function compute_surface_grad_bc(soil::SoilModel,
                                     runoff_model::CoarseGridRunoff,
                                     precip_model::AbstractPrecipModel,
                                     n̂,
                                     state⁻::Vars,
                                     diff⁻::Vars,
                                     aux⁻::Vars,
                                     t::Real
                                     )

Given a runoff model and a precipitation distribution function, compute
the surface water Neumann BC. This can be a function of time, and state.
"""
function compute_surface_grad_bc(
    soil::SoilModel,
    runoff_model::CoarseGridRunoff,
    precip_model::AbstractPrecipModel,
    n̂,
    state⁻::Vars,
    diff⁻::Vars,
    aux⁻::Vars,
    t::Real,
)
    FT = eltype(state⁻)
    incident_water_flux = precip_model(t)
    Δz = runoff_model.Δz
    water = soil.water
    param_functions = soil.param_functions
    hydraulics = water.hydraulics
    ν = param_functions.porosity

    T = get_temperature(soil.heat, aux⁻, t)
    θ_i = state⁻.soil.water.θ_i
    # Ponding Dirichlet BC
    ϑ_bc = FT(ν - θ_i)
    # Value below surface
    ϑ_below = state⁻.soil.water.ϑ_l

    # Approximate derivative of hydraulic head with respect to z
    ∂h∂z =
        FT(1) +
        (
            pressure_head(hydraulics, param_functions, ϑ_bc, θ_i) -
            pressure_head(hydraulics, param_functions, ϑ_below, θ_i)
        ) / Δz

    K =
        soil.param_functions.Ksat * hydraulic_conductivity(
            water.impedance_factor,
            water.viscosity_factor,
            water.moisture_factor,
            hydraulics,
            θ_i,
            param_functions.porosity,
            T,
            ϑ_bc / ν,# when ice is present, K still measured with ν, not νeff.
        )

    i_c = n̂ * (K * ∂h∂z)
    if incident_water_flux < -norm(i_c) # More negative if both are negative,
        #ponding BC
        K∇h⁺ = i_c
    else

        K∇h⁺ = n̂ * (-FT(2) * incident_water_flux) - diff⁻.soil.water.K∇h
    end
    return K∇h⁺
end

"""
    function compute_surface_grad_bc(soil::SoilModel,
                                     runoff_model::NoRunoff,
                                     precip_model::AbstractPrecipModel,
                                     n̂,
                                     state⁻::Vars,
                                     aux⁻::Vars,
                                     diff⁻::Vars,
                                     t::Real
                                     )

Given a runoff model and a precipitation distribution function, compute
the surface water Neumann BC. This can be a function of time, and state.
"""
function compute_surface_grad_bc(
    soil::SoilModel,
    runoff_model::NoRunoff,
    precip_model::AbstractPrecipModel,
    n̂,
    state⁻::Vars,
    diff⁻::Vars,
    aux⁻::Vars,
    t::Real,
)
    FT = eltype(state⁻)
    incident_water_flux = precip_model(t)
    K∇h⁺ = n̂ * (-FT(2) * incident_water_flux) - diff⁻.soil.water.K∇h
    return K∇h⁺
end

end


text:
Rough volatility with Python
Jim Gatheral

For Python Quants, New York, Friday May 6, 2016

Acknowledgements
The code in this iPython notebook used to be in R. I am very grateful to Yves Hilpisch and Michael Schwed for translating my R-code to Python.

For slideshow functionality I use RISE by Damián Avila.

Outline of presentation
The time series of historical volatility

Scaling properties
The RFSV model

Pricing under rough volatility

Forecasting realized variance

The time series of variance swaps

Relating historical and implied

The time series of realized variance
Assuming an underlying variance process  vs , integrated variance  1δ∫t+δtvsds  may (in principle) be estimated arbitrarily accurately given enough price data.

In practice, market microstructure noise makes estimation harder at very high frequency.
Sophisticated estimators of integrated variance have been developed to adjust for market microstructure noise. See Gatheral and Oomen [6] (for example) for details of these.
The Oxford-Man Institute of Quantitative Finance makes historical realized variance (RV) estimates freely available at http://realized.oxford-man.ox.ac.uk. These estimates are updated daily.

Each day, for 21 different indices, all trades and quotes are used to estimate realized (or integrated) variance over the trading day from open to close.
Using daily RV estimates as proxies for instantaneous variance, we may investigate the time series properties of  vt  empirically.
First load all necessary Python libraries.

In [1]:
import warnings; warnings.simplefilter('ignore')
import datetime as dt
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
from matplotlib.mlab import stineman_interp
import pandas as pd
import pandas.io.data as web
import requests
import zipfile as zi
import StringIO as sio
from sklearn import datasets, linear_model
import scipy.special as scsp
import statsmodels.api as sm
import math
import seaborn as sns; sns.set()
%matplotlib inline
Then update and save the latest Oxford-Man data.

In [2]:
url = 'http://realized.oxford-man.ox.ac.uk/media/1366/'
url += 'oxfordmanrealizedvolatilityindices.zip'
data = requests.get(url, stream=True).content
z = zi.ZipFile(sio.StringIO(data))
z.extractall()
There are many different estimates of realized variance, all of them very similar. We will use the realized kernel estimates denoted by ".rk".

In [3]:
df = pd.read_csv('OxfordManRealizedVolatilityIndices.csv', index_col=0, header=2 )
rv1 = pd.DataFrame(index=df.index)
for col in df.columns:
    if col[-3:] == '.rk':
        rv1[col] = df[col]
rv1.index = [dt.datetime.strptime(str(date), "%Y%m%d") for date in rv1.index.values]
Let's plot SPX realized variance.

In [4]:
spx = pd.DataFrame(rv1['SPX2.rk'])
spx.plot(color='red', grid=True, title='SPX realized variance',
         figsize=(16, 9), ylim=(0,0.003));

Figure 1: Oxford-Man KRV estimates of SPX realized variance from January 2000 to the current date.

In [5]:
spx.head()
Out[5]:
SPX2.rk
2000-01-03	0.000161
2000-01-04	0.000264
2000-01-05	0.000305
2000-01-06	0.000149
2000-01-07	0.000123
In [6]:
spx.tail()
Out[6]:
SPX2.rk
2016-04-27	0.000031
2016-04-28	0.000031
2016-04-29	0.000048
2016-05-02	0.000028
2016-05-03	0.000041
We can get SPX data from Yahoo using the DataReader function:

In [7]:
SPX = web.DataReader(name = '^GSPC',data_source = 'yahoo', start='2000-01-01')
SPX = SPX['Adj Close']
SPX.plot(title='SPX',figsize=(14, 8));

The smoothness of the volatility process
For  q≥0 , we define the  q th sample moment of differences of log-volatility at a given lag  Δ .( ⟨⋅⟩  denotes the sample average):

m(q,Δ)=⟨|logσt+Δ−logσt|q⟩

For example

m(2,Δ)=⟨(logσt+Δ−logσt)2⟩

is just the sample variance of differences in log-volatility at the lag  Δ .

Scaling of  m(q,Δ)  with lag  Δ
In [8]:
spx['sqrt']= np.sqrt(spx['SPX2.rk'])
spx['log_sqrt'] = np.log(spx['sqrt'])

def del_Raw(q, x):
    return [np.mean(np.abs(spx['log_sqrt'] - spx['log_sqrt'].shift(lag)) ** q)
            for lag in x]
In [9]:
plt.figure(figsize=(8, 8))
plt.xlabel('$log(\Delta)$')
plt.ylabel('$log\  m(q.\Delta)$')
plt.ylim=(-3, -.5)

zeta_q = list()
qVec = np.array([.5, 1, 1.5, 2, 3])
x = np.arange(1, 100)
for q in qVec:
    plt.plot(np.log(x), np.log(del_Raw(q, x)), 'o')
    model = np.polyfit(np.log(x), np.log(del_Raw(q, x)), 1)
    plt.plot(np.log(x), np.log(x) * model[0] + model[1])
    zeta_q.append(model[0])

print zeta_q
[0.072526605081197543, 0.14178030350291049, 0.20760692212873791, 0.27007948205423049, 0.38534332343872962]

Figure 2:  logm(q,Δ)  as a function of  logΔ , SPX.

Monofractal scaling result
From the above log-log plot, we see that for each  q ,  m(q,Δ)∝Δζq .
How does  ζq  scale with  q ?
Scaling of  ζq  with  q
In [10]:
plt.figure(figsize=(8,8))
plt.xlabel('q')
plt.ylabel('$\zeta_{q}$')
plt.plot(qVec, zeta_q, 'or')

line = np.polyfit(qVec[:4], zeta_q[:4],1)
plt.plot(qVec, line[0] * qVec + line[1])
h_est= line[0]
print(h_est)
0.131697049909

Figure 3: Scaling of  ζq  with  q .

We find the monofractal scaling relationship

ζq=qH

with  H≈0.13 .

Note however that  H  does vary over time, in a narrow range.
Note also that our estimate of  H  is biased high because we proxied instantaneous variance  vt  with its average over each day  1T∫T0vtdt , where  T  is one trading day.
Estimated  H  for all indices
We now repeat this analysis for all 21 indices in the Oxford-Man dataset.

In [11]:
def dlsig2(sic, x, pr=False):
    if pr:
        a= np.array([(sig-sig.shift(lag)).dropna() for lag in x])
        a=a ** 2
        print a.info()
    return [np.mean((sig-sig.shift(lag)).dropna() ** 2) for lag in x]
In [12]:
h = list()
nu = list()

for col in rv1.columns:
    sig = rv1[col]
    sig = np.log(np.sqrt(sig))
    sig = sig.dropna()
    model = np.polyfit(np.log(x), np.log(dlsig2(sig, x)), 1)
    nu.append(np.sqrt(np.exp(model[1])))
    h.append(model[0]/2.)

OxfordH = pd.DataFrame({'names':rv1.columns, 'h_est': h, 'nu_est': nu})

In [13]:
OxfordH
Out[13]:
h_est	names	nu_est
0	0.133954	SPX2.rk	0.321337
1	0.142315	FTSE2.rk	0.270677
2	0.113366	N2252.rk	0.320396
3	0.150251	GDAXI2.rk	0.274873
4	NaN	RUT2.rk	NaN
5	0.083370	AORD2.rk	0.359025
6	0.131013	DJI2.rk	0.317327
7	NaN	IXIC2.rk	NaN
8	0.130485	FCHI2.rk	0.291967
9	0.103370	HSI2.rk	0.281453
10	0.125847	KS11.rk	0.274713
11	0.145671	AEX.rk	0.290187
12	0.178427	SSMI.rk	0.223249
13	0.128117	IBEX2.rk	0.281662
14	0.110092	NSEI.rk	0.324883
15	0.092252	MXX.rk	0.324180
16	0.106595	BVSP.rk	0.312907
17	NaN	GSPTSE.rk	NaN
18	0.119857	STOXX50E.rk	0.337045
19	0.127094	FTSTI.rk	0.228910
20	0.133361	FTSEMIB.rk	0.298712
Distributions of  (logσt+Δ−logσt)  for various lags  Δ
Having established these beautiful scaling results for the moments, how do the histograms look?

In [14]:
def plotScaling(j, scaleFactor):
    col_name = rv1.columns[j]
    v = rv1[col_name]
    x = np.arange(1,101)

    def xDel(x, lag):
        return x-x.shift(lag)

    def sdl(lag):
        return (xDel(np.log(v), lag)).std()

    sd1 = (xDel(np.log(v), 1)).std()
    h = OxfordH['h_est'][j]
    f, ax = plt.subplots(2,2,sharex=False, sharey=False, figsize=(10, 10))

    for i_0 in range(0, 2):
        for i_1 in range(0, 2):
            la = scaleFactor ** (i_1*1+i_0*2)

            hist_val = xDel(np.log(v), la).dropna()
            std = hist_val.std()
            mean = hist_val.mean()

            ax[i_0][i_1].set_title('Lag = %s Days' %la)
            n, bins, patches = ax[i_0][i_1].hist(hist_val.values, bins=100,
                                   normed=1, facecolor='green',alpha=0.2)
            ax[i_0][i_1].plot(bins, mlab.normpdf(bins,mean,std), "r")
            ax[i_0][i_1].plot(bins, mlab.normpdf(bins,0,sd1 * la ** h), "b--")
            hist_val.plot(kind='density', ax=ax[i_0][i_1])

In [15]:
 plotScaling(1,5)

Figure 4: Histograms of  (logσt+Δ−logσt)  for various lags  Δ ; normal fit in red;  Δ=1  normal fit scaled by  Δ0.14  in blue.

Universality?
[Gatheral, Jaisson and Rosenbaum][5] compute daily realized variance estimates over one hour windows for DAX and Bund futures contracts, finding similar scaling relationships.
We have also checked that Gold and Crude Oil futures scale similarly.

Although the increments  (logσt+Δ−logσt)  seem to be fatter tailed than Gaussian.
A natural model of realized volatility
As noted originally by [Andersen et al.][1], distributions of differences in the log of realized volatility are close to Gaussian.

This motivates us to model  σt  as a lognormal random variable.
Moreover, the scaling property of variance of RV differences suggests the model:
(1)
logσt+Δ−logσt=ν(WHt+Δ−WHt)

where  WH  is fractional Brownian motion.

In [Gatheral, Jaisson and Rosenbaum][5], we refer to a stationary version of (1) as the RFSV (for Rough Fractional Stochastic Volatility) model.
Fractional Brownian motion (fBm)
Fractional Brownian motion (fBm)  {WHt;t∈R}  is the unique Gaussian process with mean zero and autocovariance function
E[WHtWHs]=12{|t|2H+|s|2H−|t−s|2H}

where  H∈(0,1)  is called the Hurst index or parameter.

In particular, when  H=1/2 , fBm is just Brownian motion.

If  H>1/2 , increments are positively correlated.% so the process is trending.

If  H<1/2 , increments are negatively correlated.% so the process is reverting.
Representations of fBm
There are infinitely many possible representations of fBm in terms of Brownian motion. For example, with  γ=12−H ,

Mandelbrot-Van Ness

WHt=CH{∫t−∞dWs(t−s)γ−∫0−∞dWs(−s)γ}.

The choice

CH=2HΓ(3/2−H)Γ(H+1/2)Γ(2−2H)−−−−−−−−−−−−−−−−−−−√

ensures that

E[WHtWHs]=12{t2H+s2H−|t−s|2H}.

Does simulated RSFV data look real?

Figure 8: Volatility of SPX (above) and of the RFSV model (below).

Remarks on the comparison
The simulated and actual graphs look very alike.

Persistent periods of high volatility alternate with low volatility periods.

H∼0.1  generates very rough looking sample paths (compared with  H=1/2  for Brownian motion).

Hence rough volatility.

On closer inspection, we observe fractal-type behavior.

The graph of volatility over a small time period looks like the same graph over a much longer time period.

This feature of volatility has been investigated both empirically and theoretically in, for example, [Bacry and Muzy][3] .

In particular, their Multifractal Random Walk (MRW) is related to a limiting case of the RSFV model as  H→0 .

Pricing under rough volatility
The foregoing behavior suggest the following model (see [Bayer et al.][2] for volatility under the real (or historical or physical) measure  P :

logσt=νWHt.

Let  γ=12−H . We choose the Mandelbrot-Van Ness representation of fractional Brownian motion  WH  as follows:

WHt=CH{∫t−∞dWPs(t−s)γ−∫0−∞dWPs(−s)γ}.

Then

==:logvu−logvtνCH{∫ut1(u−s)γdWPs+∫t−∞[1(u−s)γ−1(t−s)γ]dWPs}2νCH[Mt(u)+Zt(u)].

Note that  EP[Mt(u)|Ft]=0  and  Zt(u)  is  Ft -measurable.

To price options, it would seem that we would need to know  Ft , the entire history of the Brownian motion  Ws  for $s
Pricing under  P
Let

W~Pt(u):=2H−−−√∫utdWPs(u−s)γ

With  η:=2νCH/2H−−−√  we have  2νCHMt(u)=ηW~Pt(u)  so denoting the stochastic exponential by  E(⋅) , we may write

vu==vtexp{ηW~Pt(u)+2νCHZt(u)}EP[vu|Ft]E(ηW~Pt(u)).

The conditional distribution of  vu  depends on  Ft  only through the variance forecasts  EP[vu|Ft] ,
To price options, one does not need to know  Ft , the entire history of the Brownian motion  WPs  for $s
Pricing under  Q
Our model under  P  reads:

(2)
vu=EP[vu|Ft]E(ηW~Pt(u)).

Consider some general change of measure

dWPs=dWQs+λsds,

where  {λs:s>t}  has a natural interpretation as the price of volatility risk.

We may then rewrite (2) as

vu=EP[vu|Ft]E(ηW~Qt(u))exp{η2H−−−√∫utλs(u−s)γds}.

Although the conditional distribution of  vu  under  P  is lognormal, it will not be lognormal in general under  Q .

The upward sloping smile in VIX options means  λs  cannot be deterministic in this picture.
The rough Bergomi (rBergomi) model
Let's nevertheless consider the simplest change of measure

dWPs=dWQs+λ(s)ds,

where  λ(s)  is a deterministic function of  s . Then from (2), we would have

vu==EP[vu|Ft]E(ηW~Qt(u))exp{η2H−−−√∫ut1(u−s)γλ(s)ds}ξt(u)E(ηW~Qt(u))

where the forward variances  ξt(u)=EQ[vu|Ft]  are (at least in principle) tradable and observed in the market.

ξt(u)  is the product of two terms:

EP[vu|Ft]  which depends on the historical path $\{W_s, s
a term which depends on the price of risk  λ(s) .
Features of the rough Bergomi model
The rBergomi model is a non-Markovian generalization of the Bergomi model:
E[vu|Ft]≠E[vu|vt].

The rBergomi model is Markovian in the (infinite-dimensional) state vector  EQ[vu|Ft]=ξt(u) .
We have achieved our aim from Session 1 of replacing the exponential kernels in the Bergomi model with a power-law kernel.
We may therefore expect that the rBergomi model will generate a realistic term structure of ATM volatility skew.
Re-interpretation of the conventional Bergomi model
A conventional  n -factor Bergomi model is not self-consistent for an arbitrary choice of the initial forward variance curve  ξt(u) .

ξt(u)=E[vu|Ft]  should be consistent with the assumed dynamics.
Viewed from the perspective of the fractional Bergomi model however:

The initial curve  ξt(u)  reflects the history $\{W_s; s
The exponential kernels in the exponent of the conventional Bergomi model approximate more realistic power-law kernels.
The conventional two-factor Bergomi model is then justified in practice as a tractable Markovian engineering approximation to a more realistic fractional Bergomi model.
The stock price process
The observed anticorrelation between price moves and volatility moves may be modeled naturally by anticorrelating the Brownian motion  W  that drives the volatility process with the Brownian motion driving the price process.
Thus
dStSt=vt−−√dZt

with
dZt=ρdWt+1−ρ2−−−−−√dW⊥t

where  ρ  is the correlation between volatility moves and price moves.
Simulation of the rBergomi model
We simulate the rBergomi model as follows:

Construct the joint covariance matrix for the Volterra process  W~  and the Brownian motion  Z  and compute its Cholesky decomposition.
For each time, generate iid normal random vectors {and multiply them by the lower-triangular matrix obtained by the Cholesky decomposition} to get a  m×2n  matrix of paths of  W~  and  Z  with the correct joint marginals.
With these paths held in memory, we may evaluate the expectation under  Q  of any payoff of interest.
This procedure is very slow!

Speeding up the simulation is work in progress.
Guessing rBergomi model parameters
The rBergomi model has only three parameters:  H ,  η  and  ρ .
If we had a fast simulation, we could just iterate on these parameters to find the best fit to observed option prices. But we don't.
However, the model parameters  H ,  η  and  ρ  have very direct interpretations:

H  controls the decay of ATM skew  ψ(τ)  for very short expirations.

The product  ρη  sets the level of the ATM skew for longer expirations.

Keeping  ρη  constant but decreasing  ρ  (so as to make it more negative) pushes the minimum of each smile towards higher strikes.
So we can guess parameters in practice.
As we will see, even without proper calibration (i.e. just guessing parameters), rBergomi model fits to the volatility surface are amazingly good.
SPX smiles in the rBergomi model
In Figures 9 and 10, we show how well a rBergomi model simulation with guessed parameters fits the SPX option market as of February 4, 2010, a day when the ATM volatility term structure happened to be pretty flat.

rBergomi parameters were:  H=0.07 ,  η=1.9 ,  ρ=−0.9 .
Only three parameters to get a very good fit to the whole SPX volatility surface!
rBergomi fits to SPX smiles as of 04-Feb-2010

Figure 9: Red and blue points represent bid and offer SPX implied volatilities; orange smiles are from the rBergomi simulation.

Shortest dated smile as of February 4, 2010

Figure 10: Red and blue points represent bid and offer SPX implied volatilities; orange smile is from the rBergomi simulation.

ATM volatilities and skews
In Figures 11 and 12, we see just how well the rBergomi model can match empirical skews and vols. Recall also that the parameters we used are just guesses!

Term structure of ATM skew as of February 4, 2010

Figure 11: Blue points are empirical skews; the red line is from the rBergomi simulation.

Term structure of ATM vol as of February 4, 2010

Figure 12: Blue points are empirical ATM volatilities; the red line is from the rBergomi simulation.

Another date
Now we take a look at another date: August 14, 2013, two days before the last expiration date in our dataset.
Options set at the open of August 16, 2013 so only one trading day left.
Note in particular that the extreme short-dated smile is well reproduced by the rBergomi model.
There is no need to add jumps!
SPX smiles as of August 14, 2013

Figure 13: Red and blue points represent bid and offer SPX implied volatilities; orange smiles are from the rBergomi simulation.

The forecast formula
In the RFSV model (1),  logσt≈νWHt+C  for some constant  C .
[Nuzman and Poor][7] show that  WHt+Δ  is conditionally Gaussian with conditional expectation
E[WHt+Δ|Ft]=cos(Hπ)πΔH+1/2∫t−∞WHs(t−s+Δ)(t−s)H+1/2ds

and conditional variance

Var[WHt+Δ|Ft]=cΔ2H.

where
c=Γ(3/2−H)Γ(H+1/2)Γ(2−2H).

The forecast formula
Thus, we obtain

Variance forecast formula

(3)
EP[vt+Δ|Ft]=exp{EP[log(vt+Δ)|Ft]+2cν2Δ2H}

where

EP[logvt+Δ|Ft]=cos(Hπ)πΔH+1/2∫t−∞logvs(t−s+Δ)(t−s)H+1/2ds.

Implement variance forecast in Python
In [16]:
def c_tilde(h):
    return scsp.gamma(3. / 2. - h) / scsp.gamma(h + 1. / 2.) * scsp.gamma(2. - 2. * h)

def forecast_XTS(rvdata, h, date, nLags, delta, nu):
    i = np.arange(nLags)
    cf = 1./((i + 1. / 2.) ** (h + 1. / 2.) * (i + 1. / 2. + delta))
    ldata = rvdata.truncate(after=date)
    l = len(ldata)
    ldata = np.log(ldata.iloc[l - nLags:])
    ldata['cf'] = np.fliplr([cf])[0]
    # print ldata
    ldata = ldata.dropna()
    fcst = (ldata.iloc[:, 0] * ldata['cf']).sum() / sum(ldata['cf'])

    return math.exp(fcst + 2 * nu ** 2 * c_tilde(h) * delta ** (2 * h))
SPX actual vs forecast variance
In [17]:
rvdata = pd.DataFrame(rv1['SPX2.rk'])
nu  = OxfordH['nu_est'][0] # Vol of vol estimate for SPX
h = OxfordH['h_est'][0]
n = len(rvdata)
delta = 1
nLags = 500
dates = rvdata.iloc[nLags:n-delta].index
rv_predict = [forecast_XTS(rvdata, h=h, date=d, nLags=nLags,
                           delta=delta, nu=nu) for d in dates]
rv_actual = rvdata.iloc[nLags+delta:n].values
Scatter plot of delta days ahead predictions
In [18]:
plt.figure(figsize=(8, 8))
plt.plot(rv_predict, rv_actual, 'bo');

Figure 14: Actual vols vs predicted vols.

Superimpose actual and predicted vols
In [19]:
plt.figure(figsize=(11, 6))
vol_actual = np.sqrt(np.multiply(rv_actual,252))
vol_predict = np.sqrt(np.multiply(rv_predict,252))
plt.plot(vol_actual, "b")
plt.plot(vol_predict, "r");

Figure 15: Actual volatilities in blue; predicted vols in red.

Forecasting the variance swap curve
Finally, we forecast the whole variance swap curve using the variance forecasting formula (3).

In [20]:
def xi(date, tt, nu,h, tscale):  # dt=(u-t) is in units of years
    rvdata = pd.DataFrame(rv1['SPX2.rk'])
    return [ forecast_XTS(rvdata,h=h,date=date,nLags=500,delta=dt*tscale,nu=nu) for dt in tt]

nu = OxfordH["nu_est"][0]
h = OxfordH["h_est"][0]

def varSwapCurve(date, bigT, nSteps, nu, h, tscale, onFactor):
  # Make vector of fwd variances
  tt = [ float(i) * (bigT) / nSteps for i in range(nSteps+1)]
  delta_t = tt[1]
  xicurve = xi(date, tt, nu, h, tscale)
  xicurve_mid = (np.array(xicurve[0:nSteps]) + np.array(xicurve[1:nSteps+1])) / 2
  xicurve_int = np.cumsum(xicurve_mid) * delta_t
  varcurve1 = np.divide(xicurve_int, np.array(tt[1:]))
  varcurve = np.array([xicurve[0],]+list(varcurve1))
  varcurve = varcurve * onFactor * tscale #  onFactor is to compensate for overnight moves
  res = pd.DataFrame({"texp":np.array(tt), "vsQuote":np.sqrt(varcurve)})
  return(res)
In [21]:
def varSwapForecast(date,tau,nu,h,tscale,onFactor):
  vsc = varSwapCurve(date, bigT=2.5, nSteps=100, nu=nu, h=h,
                    tscale=tscale, onFactor=onFactor) # Creates the whole curve
  x = vsc['texp']
  y = vsc['vsQuote']
  res = stineman_interp(tau,x,y,None)

  return(res)

# Test the function

tau = (.25,.5,1,2)
date = dt.datetime(2008,9,8)
varSwapForecast(date,tau,nu=nu,h=h,tscale=252,onFactor=1)
Out[21]:
array([ 0.21949454,  0.21398188,  0.2117466 ,  0.21262899])
'Constructing a time series of variance swap curves
For each of 2,658 days from Jan 27, 2003 to August 31, 2013:

We compute proxy variance swaps from closing prices of SPX options sourced from OptionMetrics (www.optionmetrics.com) via WRDS.
We form the forecasts  EP[vu|Ft]  using (3) with 500 lags of SPX RV data sourced from The Oxford-Man Institute of Quantitative Finance (http://realized.oxford-man.ox.ac.uk).
We note that the actual variance swap curve is a factor (of roughly 1.4) higher than the forecast, which we may attribute to a combination of overnight movements of the index and the price of volatility risk.
Forecasts must therefore be rescaled to obtain close-to-close realized variance forecasts.
3-month forecast vs actual variance swaps

Figure 16: Actual (proxy) 3-month variance swap quotes in blue vs forecast in red (with no scaling factor).

Ratio of actual to forecast

Figure 17: The ratio between 3-month actual variance swap quotes and 3-month forecasts.

The Lehman weekend
Empirically, it seems that the variance curve is a simple scaling factor times the forecast, but that this scaling factor is time-varying.

We can think of this factor as having two multiplicative components: the overnight factor, and the price of volatility risk.
Recall that as of the close on Friday September 12, 2008, it was widely believed that Lehman Brothers would be rescued over the weekend. By Monday morning, we knew that Lehman had failed.
In Figure 18, we see that variance swap curves just before and just after the collapse of Lehman are just rescaled versions of the RFSV forecast curves.
We need variance swap estimates for 12-Sep-2008 and 15-Sep-2008
We proxy these by taking SVI fits for the two dates and computing the log-strips.

In [22]:
varSwaps12 =(
    0.2872021, 0.2754535, 0.2601864, 0.2544684, 0.2513854, 0.2515314,
    0.2508418, 0.2520099, 0.2502763, 0.2503309, 0.2580933, 0.2588361,
    0.2565093)

texp12 = (
    0.01916496, 0.04654346, 0.09582478, 0.19164956, 0.26830938, 0.29842574,
    0.51745380, 0.54483231, 0.76659822, 0.79397673, 1.26488706, 1.76317591,
    2.26146475)

varSwaps15 = (
    0.4410505, 0.3485560, 0.3083603, 0.2944378, 0.2756881, 0.2747838,
    0.2682212, 0.2679770, 0.2668113, 0.2706713, 0.2729533, 0.2689598,
    0.2733176)

texp15 = (
    0.01095140, 0.03832991, 0.08761123, 0.18343600, 0.26009582, 0.29021218,
    0.50924025, 0.53661875, 0.75838467, 0.78576318, 1.25667351, 1.75496235,
    2.25325120)
Actual vs predicted over the Lehman weekend
In [23]:
nu = OxfordH['nu_est'][0]
h = OxfordH['h_est'][0]
date1 = dt.datetime(2008, 9, 12)
date2 = dt.datetime(2008, 9, 15)

# Variance curve fV model forecasts
tau1000 = [ float(i) * 2.5 / 1000. for i in range(1,1001)]
vs1 = varSwapForecast(date1, tau1000, nu=nu,h=h, tscale=252, onFactor=1.29)
vs2 = varSwapForecast(date2, tau1000, nu=nu,h=h, tscale=252, onFactor=1.29)
In [24]:
plt.figure(figsize=(11, 6))
plt.plot(texp12, varSwaps12, "r")
plt.plot(texp15, varSwaps15, "b")
plt.plot(tau1000, vs1, "r--")
plt.plot(tau1000, vs2, "b--");

Figure 18: SPX variance swap curves as of September 12, 2008 (red) and September 15, 2008 (blue). The dashed curves are RFSV model forecasts rescaled by the 3-month ratio ( 1.29 ) as of the Friday close.

Remarks
We note that

The actual variance swaps curves are very close to the forecast curves, up to a scaling factor.
We are able to explain the change in the variance swap curve with only one extra observation: daily variance over the trading day on Monday 15-Sep-2008.
The SPX options market appears to be backward-looking in a very sophisticated way.
The Flash Crash
The so-called Flash Crash of Thursday May 6, 2010 caused intraday realized variance to be much higher than normal.
In Figure 19, we plot the actual variance swap curves as of the Wednesday and Friday market closes together with forecast curves rescaled by the 3-month ratio as of the close on Wednesday May 5 (which was  2.52 ).
We see that the actual variance curve as of the close on Friday is consistent with a forecast from the time series of realized variance that includes the anomalous price action of Thursday May 6.
Variance swap estimates
We again proxy variance swaps for 05-May-2010, 07-May-2010 and 10-May-2010 by taking SVI fits (see [Gatheral and Jacquier][4] ) for the three dates and computing the log-strips.

In [25]:
varSwaps5 = (
    0.4250369, 0.2552473, 0.2492892, 0.2564899, 0.2612677, 0.2659618, 0.2705928, 0.2761203,
    0.2828139, 0.2841165, 0.2884955, 0.2895839, 0.2927817, 0.2992602, 0.3116500)

texp5 = (
    0.002737851, 0.043805613, 0.120465435, 0.150581793, 0.197125257, 0.292950034,
    0.369609856, 0.402464066, 0.618754278, 0.654346338, 0.867898700, 0.900752909,
    1.117043121, 1.615331964, 2.631074606)

varSwaps7 = (
    0.5469727, 0.4641713, 0.3963352, 0.3888213, 0.3762354, 0.3666858, 0.3615814, 0.3627013,
    0.3563324, 0.3573946, 0.3495730, 0.3533829, 0.3521515, 0.3506186, 0.3594066)

texp7 = (
    0.01642710, 0.03832991, 0.11498973, 0.14510609, 0.19164956, 0.28747433, 0.36413415,
    0.39698836, 0.61327858, 0.64887064, 0.86242300, 0.89527721, 1.11156742, 1.60985626,
    2.62559890)

varSwaps10 = (
    0.3718439, 0.3023223, 0.2844810, 0.2869835, 0.2886912, 0.2905637, 0.2957070, 0.2960737,
    0.3005086, 0.3031188, 0.3058492, 0.3065815, 0.3072041, 0.3122905, 0.3299425)

texp10 = (
    0.008213552, 0.030116359, 0.106776181, 0.136892539, 0.183436003, 0.279260780,
    0.355920602, 0.388774812, 0.605065024, 0.640657084, 0.854209446, 0.887063655,
    1.103353867, 1.601642710, 2.617385352)
In [26]:
date1 = dt.datetime(2010, 5, 5)
date2 = dt.datetime(2010, 5, 7)

vsf5 = varSwapCurve(date1, bigT=2.5, nSteps=100, nu=nu, h=h, tscale=252, onFactor=2.52)
vsf7 = varSwapCurve(date2, bigT=2.5, nSteps=100, nu=nu, h=h, tscale=252, onFactor=2.52)
In [27]:
plt.figure(figsize=(11, 6))
plt.plot(texp5, varSwaps5, "r", label='May 5')
plt.plot(texp7, varSwaps7, "g", label='May 7')

plt.legend()
plt.xlabel("Time to maturity")
plt.ylabel("Variance swap quote")

plt.plot(vsf5['texp'], vsf5['vsQuote'], "r--")
plt.plot(vsf7['texp'], vsf7['vsQuote'], "g--");

Figure 19: SPX variance swap curves as of May 5, 2010 (red) and May 7, 2010 (green). The dashed curves are RFSV model forecasts rescaled by the 3-month ratio ( 2.52 ) as of the close on Wednesday May 5. The curve as of the close on May 7 is consistent with the forecast including the crazy moves on May 6.

The weekend after the Flash Crash
Now we plot forecast and actual variance swap curves as of the close on Friday May 7 and Monday May 10.

In [28]:
date1 = dt.datetime(2010,5,7)
date2 = dt.datetime(2010,5,10)

vsf7 = varSwapCurve(date1, bigT=2.5, nSteps=100, nu=nu, h=h, tscale=252, onFactor=2.52)
vsf10 = varSwapCurve(date2, bigT=2.5, nSteps=100, nu=nu, h=h, tscale=252, onFactor=2.52)
In [29]:
plt.figure(figsize=(11, 6))
plt.plot(texp7, varSwaps7, "g", label='May 7')
plt.plot(texp10, varSwaps10, "m", label='May 10')

plt.legend()
plt.xlabel("Time to maturity")
plt.ylabel("Variance swap quote")

plt.plot(vsf7['texp'], vsf7['vsQuote'], "g--")
plt.plot(vsf10['texp'], vsf10['vsQuote'], "m--");

Figure 20: The May 10 actual curve is inconsistent with a forecast that includes the Flash Crash.

Now let's see what happens if we exclude the Flash Crash from the time series used to generate the variance curve forecast.

In [30]:
plt.figure(figsize=(11, 6))
ax = plt.subplot(111)
rvdata_p = rvdata.drop((dt.datetime(2010, 5, 6)), axis=0)
rvdata.loc["2010-05-04":"2010-05-10"].plot(ax=ax, legend=False)
rvdata_p.loc["2010-05-04":"2010-05-10"].plot(ax=ax, legend=False);

Figure 21: rvdata_p has the May 6 realized variance datapoint eliminated (green line). Notice the crazy realized variance estimate for May 6!

We need a new variance curve forecast function that uses the new time series.

In [31]:
def xip(date, tt, nu,h, tscale):  # dt=(u-t) is in units of years
    rvdata = pd.DataFrame(rv1['SPX2.rk'])
    rvdata_p = rvdata.drop((dt.datetime(2010, 5, 6)), axis=0)
    return [ forecast_XTS(rvdata_p, h=h, date=date,nLags=500,
                          delta=delta_t * tscale, nu=nu) for delta_t in tt]

nu = OxfordH["nu_est"][0]
h = OxfordH["h_est"][0]

def varSwapCurve_p(date, bigT, nSteps, nu, h, tscale, onFactor):
  # Make vector of fwd variances
  tt = [ float(i) * (bigT) / nSteps for i in range(nSteps+1)]
  delta_t = tt[1]
  xicurve = xip(date, tt, nu, h, tscale)
  xicurve_mid = (np.array(xicurve[0:nSteps]) + np.array(xicurve[1:nSteps + 1])) / 2
  xicurve_int = np.cumsum(xicurve_mid) * delta_t
  varcurve1 = np.divide(xicurve_int, np.array(tt[1:]))
  varcurve = np.array([xicurve[0],]+list(varcurve1))
  varcurve = varcurve * onFactor * tscale #  onFactor is to compensate for overnight moves
  res = pd.DataFrame({"texp":np.array(tt), "vsQuote":np.sqrt(varcurve)})
  return(res)

def varSwapForecast_p(date, tau, nu, h, tscale, onFactor):
  vsc = varSwapCurve_p(date, bigT=2.5, nSteps=100, nu=nu, h=h,
                    tscale=tscale, onFactor=onFactor) # Creates the whole curve
  x = vsc['texp']
  y = vsc['vsQuote']
  res = stineman_interp(tau, x, y, None)

  return(res)

# Test the function

tau = (.25, .5 ,1, 2)
date = dt.datetime(2010, 5, 10)
varSwapForecast_p(date, tau, nu=nu, h=h, tscale=252, onFactor=1. / (1 - .35))
Out[31]:
array([ 0.26077084,  0.25255902,  0.25299844,  0.26116175])
Finally, we compare our new forecast curves with the actuals.

In [32]:
date1 = dt.datetime(2010, 5, 7)
date2 = dt.datetime(2010, 5, 10)

vsf7 = varSwapCurve(date1, bigT=2.5, nSteps=100, nu=nu, h=h, tscale=252, onFactor=2.52)
vsf10p = varSwapCurve_p(date2, bigT=2.5, nSteps=100, nu=nu, h=h, tscale=252, onFactor=2.52)
In [33]:
plt.figure(figsize=(11, 6))
plt.plot(texp7, varSwaps7, "g", label='May 7')
plt.plot(texp10, varSwaps10, "m", label='May 10')

plt.legend()
plt.xlabel("Time to maturity")
plt.ylabel("Variance swap quote")

plt.plot(vsf7['texp'], vsf7['vsQuote'], "g--")
plt.plot(vsf10p['texp'], vsf10p['vsQuote'], "m--");

Figure 22: The May 10 actual curve is consistent with a forecast that excludes the Flash Crash.

Resetting of expectations over the weekend
In Figures 20 and 22, we see that the actual variance swap curve on Monday, May 10 is consistent with a forecast that excludes the Flash Crash.
Volatility traders realized that the Flash Crash should not influence future realized variance projections.
Summary
We uncovered a remarkable monofractal scaling relationship in historical volatility.

A corollary is that volatility is not a long memory process, as widely believed.
This leads to a natural non-Markovian stochastic volatility model under  P .
The simplest specification of  dQdP  gives a non-Markovian generalization of the Bergomi model.
The history of the Brownian motion $\lbrace W_s, s
This model fits the observed volatility surface surprisingly well with very few parameters.
For perhaps the first time, we have a simple consistent model of historical and implied volatility.
References

^ Torben G Andersen, Tim Bollerslev, Francis X Diebold, and Heiko Ebens, The distribution of realized stock return volatility, *Journal of Financial Economics* **61**(1) 43-76 (2001).
^ Christian Bayer, Peter Friz and Jim Gatheral, Pricing under rough volatility, *Quantitative Finance* forthcoming, available at http://papers.ssrn.com/sol3/papers.cfm?abstract_id=2554754, (2015).
^ Emmanuel Bacry and Jean-François Muzy, Log-infinitely divisible multifractal processes, *Communications in Mathematical Physics* **236**(3) 449-475 (2003).
^ Jim Gatheral and Antoine Jacquier, Arbitrage-free SVI volatility surfaces, *Quantitative Finance* **14**(1) 59-71 (2014).
^ Jim Gatheral, Thibault Jaisson and Mathieu Rosenbaum, Volatility is rough, available at http://papers.ssrn.com/sol3/papers.cfm?abstract_id=2509457, (2014).
^ Jim Gatheral and Roel Oomen, Zero-intelligence realized variance estimation, *Finance and Stochastics* **14**(2) 249-283 (2010).
^ Carl J. Nuzman and H. Vincent Poor, Linear estimation of self-similar processes via Lamperti’s transformation, *Journal of Applied Probability* **37**(2) 429-452 (2000).


text:
using Oceananigans.Operators: Δx, Δy, ΔzC, div_xyᶜᶜᵃ, div_xzᶜᵃᶜ, div_yzᵃᶜᶜ

#####
##### Outer functions for setting velocity on boundary and filling halo beyond boundary.
#####

function fill_west_halo!(u, bc::NPBC, arch, grid, clock, state)
    @views @. u.parent[1 + grid.Hx, :, :] = 0 # fix velocity on boundary
    @launch(device(arch), config=launch_config(grid, :yz),
            _fill_west_halo!(u, bc, grid, state.velocities.v, state.velocities.w))
    return nothing
end

function fill_south_halo!(v, bc::NPBC, arch, grid, clock, state)
    @views @. v.parent[:, 1 + grid.Hy, :] = 0 # fix velocity on boundary
    @launch(device(arch), config=launch_config(grid, :xz),
            _fill_south_halo!(v, bc, grid, state.velocities.u, state.velocities.w))
    return nothing
end

function fill_bottom_halo!(w, bc::NPBC, arch, grid, clock, state)
    @views @. w.parent[:, :, 1 + grid.Hz] = 0 # fix velocity on boundary
    @launch(device(arch), config=launch_config(grid, :xy),
            _fill_bottom_halo!(w, bc, grid, state.velocities.u, state.velocities.v))
    return nothing
end

function fill_east_halo!(u, bc::NPBC, arch, grid, clock, state)
    @views @. u.parent[grid.Nx + 1 + grid.Hx, :, :] = 0 # fix velocity on boundary
    @launch(device(arch), config=launch_config(grid, :yz),
            _fill_east_halo!(u, bc, grid, state.velocities.v, state.velocities.w))
    return nothing
end

function fill_north_halo!(v, bc::NPBC, arch, grid, clock, state)
    @views @. v.parent[:, grid.Ny + 1 + grid.Hy, :] = 0 # fix velocity on boundary
    @launch(device(arch), config=launch_config(grid, :xz),
            _fill_north_halo!(v, bc, grid, state.velocities.u, state.velocities.v))
    return nothing
end

function fill_top_halo!(w, bc::NPBC, arch, grid, clock, state)
    @views @. w.parent[:, :, grid.Nz + 1 + grid.Hz] = 0 # fix velocity on boundary
    @launch(device(arch), config=launch_config(grid, :xy),
            _fill_bottom_halo!(w, bc, grid, state.velocities.u, state.velocities.v))
    return nothing
end

#####
##### Kernels that use the continuity equation to fill interior points
#####

# No penetration for x-boundaries: fill in x, loop over yz

function _fill_west_halo!(u, ::NPBC, grid, v, w)
    @loop_yz j k grid begin
        @unroll for i in 0 : -1 : (1 - grid.Hx) # integrate westward
            @inbounds u[i, j, k] = u[i+1, j, k] + Δx(i, j, k, grid) * div_yzᵃᶜᶜ(i, j, k, grid, v, w)
        end
    end
    return nothing
end

function _fill_east_halo!(u, ::NPBC, grid, v, w)
    @loop_yz j k grid begin
        @unroll for i in (grid.Nx + 2) : (grid.Nx + 1 + grid.Hx) # integrate eastward
            @inbounds u[i, j, k] = u[i-1, j, k] - Δx(i-1, j, k, grid) * div_yzᵃᶜᶜ(i-1, j, k, grid, v, w)
        end
    end
    return nothing
end

# No penetration for y-boundaries: fill in y, loop over xz

function _fill_south_halo!(v, ::NPBC, grid, u, w)
    @loop_xz i k grid begin
        @unroll for j in 0 : -1 : (1 - grid.Hy) # integrate southward
            @inbounds v[i, j, k] = v[i, j+1, k] + Δy(i, j, k, grid) * div_xzᶜᵃᶜ(i, j, k, grid, u, w)
        end
    end
    return nothing
end

function _fill_north_halo!(v, ::NPBC, grid, u, w)
    @loop_xz i k grid begin
        @unroll for j in (grid.Ny + 2) : (grid.Ny + 1 + grid.Hy) # integrate northward
            @inbounds v[i, j, k] = v[i, j-1, k] - Δy(i, j-1, k, grid) * div_xzᶜᵃᶜ(i, j-1, k, grid, u, w)
        end
    end
    return nothing
end

# No penetration for z-boundaries: fill in z, loop over xy

function _fill_bottom_halo!(w, ::NPBC, grid, u, v)
    @loop_xy i j grid begin
        @unroll for k in 0 : -1 : (1 - grid.Hz) # integrate downwards
            @inbounds w[i, j, k] = w[i, j, k+1] + ΔzC(i, j, k, grid) * div_xyᶜᶜᵃ(i, j, k, grid, u, v)
        end
    end
    return nothing
end

function _fill_top_halo!(w, ::NPBC, grid, u, v)
    @loop_xy i j grid begin
        @unroll for k in (grid.Nz + 2) : (grid.Nz + 1 + grid.Hz) # integrate upwards
            @inbounds w[i, j, k] = w[i, j, k-1] - ΔzC(i, j, k-1, grid) * div_xyᶜᶜᵃ(i, j, k-1, grid, u, v)
        end
    end
    return nothing
end


text:
using Test
using LogicCircuits

include("helper/plain_logic_circuits.jl")

@testset "Smooth test" begin
    for file in [zoo_sdd_file("random.sdd")] #  save some test time;zoo_psdd_file("plants.psdd"),
        c1 = load_logic_circuit(file)
        c2 = smooth(c1)
        @test model_count(c1) == model_count(c2)
        @test !issmooth(c1)
        @test issmooth(c2)
        @test c1 !== c2
        @test smooth(c2) === c2
    end
end

@testset "Forget test" begin
    for file in [zoo_sdd_file("random.sdd")] #  save some test time;zoo_psdd_file("plants.psdd"),
        c1 = load_logic_circuit(file)
        vars = variables(c1)
        c2 = forget(c1, x -> x > maximum(vars))
        @test c2 === c1
        c3 = forget(c1, x -> x > 10)
        @test c3 !== c2
        @test variables(c3) == BitSet(1:10)
        c4 = propagate_constants(c1)
        c5 = forget(c4, x -> x == 10)
        @test !(10 in variables(c5))
    end

end

@testset "Propagate constants test" begin
    for file in [zoo_sdd_file("random.sdd")] #  save some test time;zoo_psdd_file("plants.psdd"),
        c1 = load_logic_circuit(file)
        c2 = propagate_constants(c1)
        @test model_count(c1) == model_count(c2)
        (false_node, true_node) = canonical_constants(c2)
        @test false_node === nothing
        @test true_node === nothing
        c3 = propagate_constants(c2)
        @test c3 === c2
    end
end

@testset "Deepcopy test" begin
    n0 = little_3var()
    n1 = deepcopy(n0, 0)
    @test n0 === n1
    n2 = deepcopy(n0, 1)
    @test n2 !== n0
    @test all(children(n2) .=== children(n0))
    n3 = deepcopy(n0, 2)
    @test n3 !== n0
    @test all(children(n3) .!== children(n0))
    @test all(vcat(children.(children(n3))...) .=== vcat(children.(children(n0))...))

    for depth in [5, typemax(Int)]
        n3 = load_logic_circuit(zoo_sdd_file("random.sdd"))
        n4 = deepcopy(n3, depth)
        @test num_nodes(n3) == num_nodes(n4)
        @test num_edges(n3) == num_edges(n4)
        @test model_count(n3) == model_count(n4)
    end
end

@testset "Condition test" begin
    c1 = load_logic_circuit(zoo_sdd_file("random.sdd"))

    lit = Lit(num_variables(c1) + 1)
    @test c1 === condition(c1, lit)

    lit1 = Lit(1)
    c2 = condition(c1, lit1)
    dict = canonical_literals(c2)
    @test haskey(dict, lit1)
    @test !haskey(dict, -lit1)
    c3 = condition(c2, -lit1)
    @test isfalse(c3)
    c4 = condition(c2, lit1)
    @test c4 == c2

    c1 = little_2var()
    c2 = condition(c1, lit1)
    @test num_nodes(c2) == 6
    @test num_edges(c2) == 5

    lit2 = Lit(2)
    c3 = condition(c1, lit2)
    @test num_nodes(c3) == 6
    @test num_edges(c3) == 6

    c4 = condition(c2, lit2)
    @test num_nodes(c4) == 4
    @test num_edges(c4) == 3

    lit1 = compile(PlainLogicCircuit, Lit(1))
    litn1 = compile(PlainLogicCircuit, - Lit(1))
    lit2 = compile(PlainLogicCircuit, Lit(2))
    or = Plain⋁Node([lit1])
    c5 = or & lit2 & litn1
    @test isfalse(condition(c5, - Lit(1)))
end

@testset "Split test" begin
    or = c0 = little_5var()
    and = children(or)[1]
    v = Var(1)
    c1, _ = split(c0, (or, and), v)
    @test num_nodes(c1) == 22
    @test num_edges(c1) == 24
    @test c1.children[1].children[2] === c1.children[2].children[2] === c0.children[1].children[2]
    @test c1.children[1].children[1].children[2] === c1.children[2].children[1].children[2] === c0.children[1].children[1].children[2]

    or = c1.children[1].children[1]
    and = children(or)[2]
    v = Var(4)
    c2, _ = split(c1, (or, and), v)
    @test num_nodes(c2) == 24
    @test num_edges(c2) == 29
    n1 = c2.children[2].children[1].children[2]
    n2 = c2.children[1].children[1].children[2]
    @test n1 in c1
    @test !(n2 in c1)
    @test n1.children[1] == n2.children[1]

    or = c0 = little_5var()
    and = children(c0)[1]
    c1, _ = split(c0, (or, and), v; depth=0)
    c2, _ = split(c0, (or, and), v; depth=1)
    c3, _ = split(c0, (or, and), v; depth=2)
    c4, _ = split(c0, (or, and), v; depth=3)
    @test num_nodes(c1) == num_nodes(c2) < num_nodes(c3) < num_nodes(c4)
end

@testset "Merge test" begin
    n0 = little_3var()
    or1 = n0.children[1].children[2]
    or2 = n0.children[2].children[2]
    n1 = merge(n0, or1, or2)
    @test n0.children[1].children[2] != n0.children[2].children[2]
    @test n1.children[1].children[2] == n1.children[2].children[2]
    @test num_nodes(n1) == (num_nodes(n0) - 1)
end

@testset "Clone test" begin
    lit1 = compile(PlainLogicCircuit, Lit(1))
    litn1 = compile(PlainLogicCircuit, - Lit(1))
    lit2 = compile(PlainLogicCircuit, Lit(2))
    litn2 = compile(PlainLogicCircuit, - Lit(2))
    or = lit1 | litn1
    and1 = lit2 & or
    and2 = or & litn2
    c1 = root = and1 | and2
    c2 = clone(c1, and1, and2, or)
    @test num_nodes(c2) == num_nodes(c1) + 1
    @test num_edges(c2) == num_edges(c1) + 2
    @test c2.children[1] === and1
    @test c2.children[2] !== and2
    @test c2.children[2].children[2] === and2.children[2]
    @test c2.children[2].children[1] !== or
    @test c2.children[1].children[2] === or
    @test all(c2.children[2].children[1].children .=== or.children)
end

@testset "Random splits" begin
    c0 = smooth(load_logic_circuit(zoo_sdd_file("random.sdd")))
    @test issmooth(c0)
    @test isdecomposable(c0)
    for _ in 1 : 4
        for depth in 0 : 3
            c1, _ = split_step(c0; loss=random_split, depth=depth, sanity_check=true)
            @test issmooth(c1)
            @test isdecomposable(c1)
            c0 = c1
            # TODO
            # @test isstructdecomposable(c1)
            # @test isdeterministic(c1)
        end
    end
    c1 = struct_learn(c0; maxiter=10, stop=(x-> num_edges(x) > num_edges(c0) + 10))
    @test issmooth(c1)
    @test isdecomposable(c1)
    # TODO
    # @test isstructdecomposable(c1)
    # @test isdeterministic(c1)
end

@testset "Clone Candidates" begin
    lit1 = compile(PlainLogicCircuit, Lit(1))
    litn1 = compile(PlainLogicCircuit, - Lit(1))
    lit2 = compile(PlainLogicCircuit, Lit(2))
    litn2 = compile(PlainLogicCircuit, - Lit(2))
    or = lit1 | litn1
    and1 = lit2 & or
    and2 = or & litn2
    c1 = root = and1 | and2

    cands = clone_candidates(c1)
    candidates = [(k, v[1], v[2]) for (k,v) in cands]

    @assert (or, and1, and2) in candidates
end

@testset "Split regression test 1" begin
    vtree = Vtree(784, :balanced)
    c0 = fully_factorized_circuit(StructLogicCircuit,vtree)
    c1, _ = split_step(c0; loss=random_split, depth=2, sanity_check=true)
    @test issmooth(c1)
    @test isdecomposable(c1)
end

@testset "Split candidate test" begin
    a = compile(PlainLogicCircuit, Lit(1))
    b = compile(PlainLogicCircuit, Lit(2))
    c = compile(PlainLogicCircuit, Lit(3))
    d = compile(PlainLogicCircuit, Lit(4))
    e = compile(PlainLogicCircuit, -Lit(4))

    f = disjoin(conjoin(a, b), conjoin(c, d))
    candidates, variable_scope = split_candidates(f)
    @test length(candidates) == 0

    g = disjoin(conjoin(disjoin(conjoin(d, b), conjoin(e, b)), a))
    candidates, variable_scope = split_candidates(g)
    @test length(candidates) == 1
    @test candidates[1][1] == g
end

@testset "Circuit standardize test" begin
    a = compile(PlainLogicCircuit, Lit(1))
    b = compile(PlainLogicCircuit, Lit(2))
    c = compile(PlainLogicCircuit, Lit(3))
    d = compile(PlainLogicCircuit, Lit(4))

    e = conjoin(a, b, c)
    circuit = standardize_circuit(e)
    @test circuit.children[1].literal == Lit(1)
    @test circuit.children[2].children[1].children[1].literal == Lit(2)
    @test circuit.children[2].children[1].children[2].literal == Lit(3)

    e = disjoin(a, b)
    f = disjoin(c, d)
    g = disjoin(e, f)
    circuit = standardize_circuit(g)
    @test circuit.children[1].literal == Lit(1)
    @test circuit.children[2].literal == Lit(2)
    @test circuit.children[3].literal == Lit(3)
    @test circuit.children[4].literal == Lit(4)
end


text:
################################################################################
#
#  Naive relation search: Based on coefficient size only
#
################################################################################

function class_group_random_ideal_relation(clg::ClassGrpCtx, r::Int,
                                           I::NfOrdIdl = rand(clg.FB.ideals))
  s = 1
  if r < 2
    r = 2
  end
  for i = 1:r
    I = I*rand(clg.FB.ideals)
    I, g = reduce_ideal_class(I)
    s *= g
  end
  return s;
end

# Do better: re-use partial hnf, check rank mod p, ...

################################################################################
#
#  Relation matrix processing
#
################################################################################

# a large protion is now outsourced to LinearAlgebra/Module.jl

function class_group_process_relmatrix(clg::ClassGrpCtx)

  @vtime_add_elapsed :ClassGroup 1 clg :hnf_time clg.h = check_index(clg.M)
  clg.hnf_call += 1
end

function class_group_get_pivot_info(clg::ClassGrpCtx)
  # Extract information about (missing) pivot elements.
  # If we are in the full rank case, they come from the hnf itself,
  # Otherwise we look at the echelon form of the reduction.

  @vtime_add_elapsed :ClassGroup 1 clg :hnf_time h = check_index(clg.M)
  clg.h = h
  return (h, non_trivial_pivot(clg.M)::IntSet)
end

################################################################################
#
#  Main loop to find relations
#
################################################################################

function class_group_find_relations(clg::ClassGrpCtx; val = 0, prec::Int = 100,
                limit::Int = 10)
  clg.hnf_call = 0
  clg.rel_cnt = 0
  clg.bad_rel = 0

  n = degree(clg.FB.ideals[1].parent.order)
  t = time_ns()
  I = []
  O = parent(clg.FB.ideals[1]).order
  sqrt_disc = isqrt(abs(discriminant(O)))
  sqrt_disc = max(sqrt_disc, 1000)
  np = nbits(sqrt_disc)+30

  f = 0

  for i in clg.FB.ideals
    f = class_group_small_real_elements_relation_start(clg, i, limit = limit,
                        prec = prec, val = val)

    push!(I, f)
    f.vl = val
    while true
      e = class_group_small_real_elements_relation_next(I[end])
      n = abs(norm_div(e, norm(I[end].A), np))
#      if n==0 || e==0
##        println("found ", e, " of norm ", n)
#      end
#        print_with_color(:blue, "norm OK:")
#        println(n//norm(I[end].A), " should be ", sqrt_disc)
      if nbits(numerator(n)) > np-10
#        prec = Int(ceil(prec*1.2))
#        print_with_color(:red, "norm too large:")
#        println(n, " should be ", sqrt_disc)
#        println("offending element is ", e)
#        println("offending ideal is ", I[end].A)
#        println("skipping ideal (for now)")
        break
      end
      f = class_group_add_relation(clg, e, n, norm(I[end].A), integral = true)
      if f
        I[end].cnt += 1
        break
      else
        I[end].bad += 1
        if I[end].bad > (clg.bad_rel/clg.rel_cnt)*2
          @v_do :ClassGroup 2 println("too slow in getting s.th. for ", i,
                          "\ngood: ", I[end].cnt,  " bad: ",  I[end].bad,
                          " ratio: ", (clg.bad_rel/clg.rel_cnt))
          break
        end
      end
    end
    @v_do :ClassGroup_gc 1 gc()
    #want 10% more relations than ideals...
    if cols(clg.M) *1.1 < rows(clg.M)
      @vprint :ClassGroup 1 println("rel mat probably full rank, leaving phase 1");
      while length(I) < length(clg.FB.ideals)
        push!(I, class_group_small_real_elements_relation_start(clg, clg.FB.ideals[length(I)+1], limit = limit, prec = prec, val = val))
      end
      break
    end
  end

  @v_do :ClassGroup 1 println("used ", (time_ns()-t)/1e9,
                  " sec for small elts, so far ", clg.time[:hnf_time]/1e9,
                  " sec for hnf in ", clg.hnf_call, " calls");
  @v_do :ClassGroup 1 println("added ", clg.rel_cnt, " good relations and ",
                  clg.bad_rel, " bad ones, ratio ", clg.bad_rel/clg.rel_cnt)

  @v_do :ClassGroup 1 println("total ", rows(clg.M), " relations (incl. orbit)")

  s = time_ns()

  class_group_process_relmatrix(clg)
  h, piv = class_group_get_pivot_info(clg)

  println("current h=$h")

  idl = clg.FB.ideals
  want_extra = 5
  bad_h = false
  bad_norm = 0
  new_rel = 0
  while h != 1 && (h==0 || want_extra > 0)
    for i in sort([ x for x in piv], lt = >)
      E = I[i]
      lt = max(100, round((clg.bad_rel/clg.rel_cnt)*2))

      while true
        if (E.cnt==0 && E.bad > lt) || (E.cnt != 0 && (bad_h || E.bad > lt))
          @v_do :ClassGroup 2 println("cnt ", E.cnt, " bad ", E.bad, " limit ", lt)
          @v_do :ClassGroup 2 println("re-starting (at ", i, ") ")
          @v_do :ClassGroup 3 println("for ideal ", E.A)

          A = idl[i]
          while norm(A) < sqrt_disc
            A *= rand(idl)
          end
          bad_norm = 0

          I[i] = class_group_small_real_elements_relation_start(clg, A,
                          val = E.vl, limit = limit, prec = prec)
          E = I[i]
        end
        e = class_group_small_real_elements_relation_next(E)
        n = abs(norm_div(e, norm(E.A), np))

        if nbits(numerator(n)) > np-10
          bad_norm += 1
          if bad_norm /(E.cnt + E.bad + 1) > 0.1
#            print_with_color(:red, "too many large norms, changing ideal\n")
#            println("bad_norm: ", bad_norm, " cnt: ", E.cnt, " bad ", E.bad)
            A = idl[i]
            while norm(A) < sqrt_disc
              A *= rand(idl)
            end
            I[i] = class_group_small_real_elements_relation_start(clg, A,
                            val = E.vl, limit = limit, prec = prec)
            E = I[i]
            bad_norm = 0
          end
          #= CF: think careful here
           - norm might be wrong as we did not use enough primes
           - use as large prime variant
           - bad chance for smooth
           lets skip it for now
          =#
          continue;
        end
        if class_group_add_relation(clg, e, n, norm(E.A), orbit = false, integral = true)
          E.cnt += 1
          new_rel += 1
          break
        else
          E.bad += 1
        end
        if  clg.bad_rel - clg.last > 1000000
#          global AllRels = (i, I[i], E)
          error("to bad in finding rel")
        end
      end
      if new_rel > 20
        break
      end
    end
    last_h = h
    l_piv = piv
    new_rel = 0
    class_group_process_relmatrix(clg)
    h, piv = class_group_get_pivot_info(clg)
    println("current h=$h from $(clg.M)")

    if h != 0
      if h==1
        return h, piv
      end
      @v_do :ClassGroup 1 println("full rank: current h = ", h,
                      " want ", want_extra, " more")
      if h == last_h
        want_extra -= 1
      else
        want_extra = 15
      end
    end
    if length(l_piv) - length(piv) < length(l_piv)/2
      bad_h = true
    else
      bad_h = false
    end
    @v_do :ClassGroup_gc 1 gc()
  end

  @v_do :ClassGroup 1 println("used ", (time_ns()-s)/1e9, " total so far ",
                  clg.time[:hnf_time]/1e9, " sec for hnf in ", clg.hnf_call, " calls");
  @v_do :ClassGroup 1 println("added ", clg.rel_cnt, " good relations and ",
                  clg.bad_rel, " bad ones, ratio ", clg.bad_rel/clg.rel_cnt)

  class_group_process_relmatrix(clg)
  h, piv = class_group_get_pivot_info(clg)

  return h, piv
end

function class_group_find_relations2(clg::ClassGrpCtx; val = 0, prec = 100,
                limit = 10)
  clg.hnf_call = 0
  clg.rel_cnt = 0
  clg.bad_rel = 0

  n = degree(clg.FB.ideals[1].parent.order)
  t = time_ns()
  I = []
  O = parent(clg.FB.ideals[1]).order
  sqrt_disc = isqrt(abs(discriminant(O)))
  sqrt_disc = max(sqrt_disc, 1000)
  np = nbits(sqrt_disc)+30

  local f

  old_r = rank(clg.M)
#  println("starting with $old_r in rank increase")

  nI = length(clg.FB.ideals)
  Idl = clg.FB.ideals
  for i in nI:-1:1
    I = Idl[i]
    too_slow = false
    f = class_group_small_real_elements_relation_start(clg, I,
                                       limit = limit, prec = prec, val = val)

    f.vl = val
    while true
      e = class_group_small_real_elements_relation_next(f)
      n = abs(norm_div(e, norm(f.A), np))
      if nbits(numerator(n)) > np-10 || f.restart > 0
#        print_with_color(:red, "norm too large or restarting: $(f.restart)")
#        println(n, " should be ", sqrt_disc)
#        println("offending element is ", e)
#        println("skipping ideal (for now)")
        break
      end
      fl = class_group_add_relation(clg, e, n, norm(f.A), orbit = false)
      if fl
        f.cnt += 1
        if f.cnt % 20 == 0
          a = rank(clg.M)
          if (a-old_r) < 0.5
            @v_do :ClassGroup 2 println("rank too slow $a ($old_r) and $(clg.rel_mat_full_rank)")
            too_slow = true
            break
          end
          old_r = a
        else
          break
        end
      else
        f.bad += 1
        if clg.bad_rel > clg.rel_cnt && f.bad > max(10, (clg.bad_rel/clg.rel_cnt)*2)
          @v_do :ClassGroup 2 println("too slow in getting s.th. for ", i,
                          "\ngood: ", f.cnt,  " bad: ",  f.bad,
                          " ratio: ", (clg.bad_rel/clg.rel_cnt))
          too_slow = true
          break
        end
      end
    end
    @v_do :ClassGroup_gc 1 gc()
    if too_slow
      break
    end
  end

  @v_do :ClassGroup 1 println("used ", (time_ns()-t)/1e9,
                  " sec for small elts, so far ", clg.time[:hnf_time]/1e9,
                  " sec for hnf in ", clg.hnf_call, " calls");
  @v_do :ClassGroup 1 println("added ", clg.rel_cnt, " good relations and ",
                  clg.bad_rel, " bad ones, ratio ", clg.bad_rel/clg.rel_cnt)

  s = time_ns()

  h, piv = class_group_get_pivot_info(clg)

  @vprint :ClassGroup 1 "Target rank: $(length(clg.FB.ideals))\nCurrent rank: $(rank(clg.M))\nTentative class number: $(h)"

  want_extra = 5
  bad_h = false
  no_rand = 1
  a_old = rank(clg.M)
  while h != 1 && (h==0 || want_extra > 0)
    for i in sort([ x for x in piv], lt = >)
      I = Idl[i]
      lt = max(100, round((clg.bad_rel/clg.rel_cnt)*2))

      while true
        A = Idl[i]
        j = 0
        while j < no_rand #  && norm(A) < sqrt_disc
          A *= rand(Idl)
          j += 1
        end
        bad_norm = 0
#        println("using ideal of norm $(norm(A)) no_rand $no_rand")

        E = 0
        E = class_group_small_real_elements_relation_start(clg, A,
                        val = val, limit = limit, prec = prec)
        no_rand_local = no_rand
        while true
          e = class_group_small_real_elements_relation_next(E)
          n = abs(norm_div(e, norm(E.A), np))
          if nbits(numerator(n)) > np-10 || E.restart > 5
#            @v_do :ClassGroup 2 begin
#              print_with_color(:red, "2:norm too large (or restarting):")
#              println(n, " should be ", sqrt_disc)
#              println("offending element is ", e)
#              println("prec now ", prec)
#            end
            A = Idl[i]
            j = 0
            # TH: without added no_rand_local < nI it crashes sometimes
            #     but I don't know what I am doing
            while norm(A) < sqrt_disc && j < no_rand_local && no_rand_local < nI
              A *= rand(Idl[(nI-no_rand_local):nI])
              j += 1
            end
            no_rand_local = min(nI-1, no_rand_local+1)
#            println("using ideal $A of norm $(norm(A))")
            E = class_group_small_real_elements_relation_start(clg, A,
                            val = E.vl, limit = limit, prec = prec)
            #= CF: think careful here
             - norm might be wrong as we did not use enough primes
             - use as large prime variant
             - bad chance for smooth
             lets skip it for now
            =#
            continue;
          end
          if class_group_add_relation(clg, e, n, norm(E.A))
            E.cnt += 1
            #print_with_color(:green, "success\n")
            break
            if length(clg.RS) % 20 == 0
              #print_with_color(:blue, "found rels, trying again\n")
            end
          else
            E.bad += 1
          end
        end
        if length(clg.RS) % 20 == 0 # cols(clg.M)*0.1
#          print_with_color(:blue, "2:found rels, trying again\n")
          break
        end
      end
      if length(clg.RS) % 20  == 0# cols(clg.M)*0.1
#        print_with_color(:blue, "3:found rels, trying again\n")
        break
      end
    end
    last_h = h
    l_piv = piv
    last_rank = rank(clg.M)
    class_group_process_relmatrix(clg)
    a = rank(clg.M)
#    println("rank increase gave $a was $a_old")
    h, piv = class_group_get_pivot_info(clg)
    if (a-a_old) < 0.5 * clg.bad_rel/(length(clg.RS) + clg.bad_rel)
      no_rand += 5
      no_rand = min(no_rand, length(clg.FB.ideals))
      no_rand = min(no_rand, Int(floor(length(clg.FB.ideals)*0.1)))
    end
    a_old = a
    if h != 0
      if h==1
        return h, piv
      end
      @vprint :ClassGroup 1 "Now have $(clg.M)"
      @v_do :ClassGroup 1 println("full rank: current h = ", h,
                      " want ", want_extra, " more")
      if h == last_h
        want_extra -= 1
      else
        want_extra = 15
      end
    end
    if length(l_piv) - length(piv) < length(l_piv)/2
      bad_h = true
    else
      bad_h = false
    end
    @v_do :ClassGroup_gc 1 gc()
  end

  @v_do :ClassGroup 1 println("used ", (time_ns()-s)/1e9, " total so far ",
                  clg.time[:hnf_time]/1e9, " sec for hnf in ", clg.hnf_call, " calls");
  @v_do :ClassGroup 1 println("added ", clg.rel_cnt, " good relations and ",
                  clg.bad_rel, " bad ones, ratio ", clg.bad_rel/clg.rel_cnt)
  class_group_process_relmatrix(clg)
  h, piv = class_group_get_pivot_info(clg)
end

# add one/ a few more relations

function class_group_find_new_relation(clg::ClassGrpCtx; val = 0, prec = 100,
                limit = 10, extra = 1)
  if !isdefined(clg, :randomClsEnv)
    clg.randomClsEnv = random_init(clg.FB.ideals)
  end

  O = parent(clg.FB.ideals[1]).order
  sqrt_disc = isqrt(abs(discriminant(O)))
  sqrt_disc = max(sqrt_disc, 1000)
  np = nbits(sqrt_disc)+30

  rat = max(clg.bad_rel/clg.rel_cnt, 20.0)

  while true
    I = random_get(clg.randomClsEnv)
#    println("trying in ideal $I");
    E = class_group_small_real_elements_relation_start(clg, I,
                            val = val, limit = limit, prec = prec)

    while true
      e = class_group_small_real_elements_relation_next(E)
      n = abs(norm_div(e, norm(E.A), np))
      if nbits(numerator(n)) > np-10 || E.restart > 2
        break;
      end
      if class_group_add_relation(clg, e, n, norm(E.A))
        E.cnt += 1
        extra -= 1
        if extra <= 0
          return
        end
      else
        E.bad += 1
        if E.bad > rat
          break
        end
      end
    end
  end
  class_group_process_relmatrix(clg)
end


text:
using SparseIR
using Test
using LinearAlgebra

import PyCall: pyimport, PyNULL, PyVector

@testset "poly.getitem" begin
    lambda = 10.0
    beta = 1.0
    basis = FiniteTempBasis(LogisticKernel(lambda), fermion, beta)
    u = basis.u
    @test u[1](0.0) ≈ u.o.__getitem__(0)(0.0)
end

@testset "poly.overlap" begin
    lambda = 10.0
    beta = 1.0
    basis = FiniteTempBasis(LogisticKernel(lambda), fermion, beta)
    f = x -> 2 * x
    @test overlap(basis.u[1], f) isa Float64

    overlap_res = Matrix{Float64}(undef, size(basis), size(basis))
    for j in 1:size(basis), i in 1:size(basis)
        overlap_res[i, j] = overlap(basis.u[i], basis.u[j])
    end
    @test overlap_res ≈ Matrix(I, size(basis), size(basis))

    # Some topy spectral function
    gaussian(x, mu, sigma) = exp.(-((x .- mu) / sigma) .^ 2) / (sqrt(π) * sigma)
    function rho(omega)
        return 0.2 * gaussian(omega, 0.0, 0.15) +
               0.4 * gaussian(omega, 1.0, 0.8) + 0.4 * gaussian(omega, -1.0, 0.8)
    end

    rhol = overlap(basis.v, rho)
end