Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add filer and smoothing for regime switching model #250

Open
wants to merge 17 commits into
base: master
Choose a base branch
from

Conversation

Shunsuke-Hori
Copy link
Collaborator

@Shunsuke-Hori Shunsuke-Hori commented Oct 8, 2019

Filtering and smoothing features for regime switching model are added. Some tests for them are added too.

Copy link
Member

@cc7768 cc7768 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This looks cool! Thanks for adding this @Shunsuke-Hori

I haven't walked through and double checked the math, but have added some code comments.

Are the tests based on other code that make you relatively confident they're correct?

@@ -85,3 +85,206 @@ function hamilton_filter(y::AbstractVector, h::Integer)
y_trend = y - y_cycle
return y_cycle, y_trend
end

"""
- `g::Function`: Parameterized `Function` that recieves three arguments
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We should break these lines a little more frequently than they're currently broken -- I don't remember what the maximum number of characters per line we decided to enforce for Jula though.

@sglyon should remember.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@cc7768 It totally slipped my mind. I think it is 80.

src/filter.jl Outdated
mutable struct RegimeSwitchingModel{TD <: AbstractArray, TP <: AbstractMatrix, TPS <: AbstractMatrix}
g::Function
y::TD
parameter
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there a reason that we don't want this typed? By having a mutable structure and not having this argument typed, we can change it to anything we want -- For example, I can set rsm.parameter = 1 and it won't complain. Is there a case where this changes types?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there a reason that we don't want this typed?

Kind of, I thought some may want to use user-defined struct like this rather than NamedTuple. But I also think it should be typed. Maybe I can force it to be NamedTuple.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Perfect -- This is a use case Julia should handle very cleanly! You can add a type parameter to RegimeSwitchingModel and then just refer to that type parameter when typing the argument. I'd make suggestion, but I can only figure out how to suggest one line change, but here's what it would look like:

mutable struct RegimeSwitchingModel{TD <: AbstractArray, TP <: AbstractMatrix, TPS <: AbstractMatrix, TPRM}
g::Function
y::TD
parameter::TPRM
...

This will allow the functions to specialize to the user's type (there will be a new type of RegimeSwitchingModel for every type that gets put in as a parameter) and still prevent people from doing "weird" changes to the parameter data.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh I see, that makes sense and is way better. Thank you!!

src/filter.jl Outdated
period 0.

##### Returns
- `p`: likelihood at each period
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It wouldn't hurt to emphasize that the returns are each of size T x M where T is the number of observations and M is the number of regimes.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's fair. I'll do so.

T = size(y, 1)
p = Vector{Float64}(undef, T)
logL = 0
for t in 1:T
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I like how the code is structured that allows you to write such a simple (and beautiful!) for-loop to do the core of the computation.

src/filter.jl Outdated
rsm.prob_smoothed[t, s] = sum(p_s_joint_smoothed[t, s, :])
end
end
verbose ? (return p_s_joint_smoothed) : (return nothing)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Does this if statement in return cause confusion to the compiler? This could prevent it from inferring the type of the output of this function. Is there a reason to not just return p_s_joint_smoothed?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think it does create a problem -- It identifies it as a union of nothing and the array type:

julia> @code_warntype test(rsm, true)
Variables
  #self#::Core.Compiler.Const(test, false)
  rsm::Main.QuantEcon.RegimeSwitchingModel{Array{Float64,2},Array{Float64,2},Array{Float64,2}}
  verbose::Bool

Body::Union{Nothing, Array{Float64,1}}
1 ─ %1 = Main.randn(10)::Array{Float64,1}
│   %2 = (:verbose,)::Core.Compiler.Const((:verbose,), false)
│   %3 = Core.apply_type(Core.NamedTuple, %2)::Core.Compiler.Const(NamedTuple{(:verbose,),T} where T<:Tuple, false)
│   %4 = Core.tuple(verbose)::Tuple{Bool}
│   %5 = (%3)(%4)::NamedTuple{(:verbose,),Tuple{Bool}}
│   %6 = Main.QuantEcon.smooth!::Core.Compiler.Const(Main.QuantEcon.smooth!, false)
│   %7 = Core.kwfunc(%6)::Core.Compiler.Const(getfield(Main.QuantEcon, Symbol("#kw##smooth!"))(), false)
│   %8 = Main.QuantEcon.smooth!::Core.Compiler.Const(Main.QuantEcon.smooth!, false)
│   %9 = (%7)(%5, %8, rsm, %1)::Union{Nothing, Array{Float64,1}}
└──      return %9

julia> @code_warntype test(rsm, false)
Variables
  #self#::Core.Compiler.Const(test, false)
  rsm::Main.QuantEcon.RegimeSwitchingModel{Array{Float64,2},Array{Float64,2},Array{Float64,2}}
  verbose::Bool

Body::Union{Nothing, Array{Float64,1}}
1 ─ %1 = Main.randn(10)::Array{Float64,1}
│   %2 = (:verbose,)::Core.Compiler.Const((:verbose,), false)
│   %3 = Core.apply_type(Core.NamedTuple, %2)::Core.Compiler.Const(NamedTuple{(:verbose,),T} where T<:Tuple, false)
│   %4 = Core.tuple(verbose)::Tuple{Bool}
│   %5 = (%3)(%4)::NamedTuple{(:verbose,),Tuple{Bool}}
│   %6 = Main.QuantEcon.smooth!::Core.Compiler.Const(Main.QuantEcon.smooth!, false)
│   %7 = Core.kwfunc(%6)::Core.Compiler.Const(getfield(Main.QuantEcon, Symbol("#kw##smooth!"))(), false)
│   %8 = Main.QuantEcon.smooth!::Core.Compiler.Const(Main.QuantEcon.smooth!, false)
│   %9 = (%7)(%5, %8, rsm, %1)::Union{Nothing, Array{Float64,1}}
└──      return %9

It does add "complexity" relative to return p_s_joint_smoothed though, so is there a reason not to return it always?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

is there a reason not to return it always?

No, when I realized p_s_joint_smoothed should be returned, I left return nothing case so that I didn't need to make change to existing parts, but that was a silly idea. Now I agree with you, it's just confusing. I'll fix it.

@Shunsuke-Hori
Copy link
Collaborator Author

@cc7768 Thank you for your review!! Your comments improve the codes a lot.

The test is based on python statsmodels, which I assumed correct. (https://www.statsmodels.org/dev/examples/notebooks/generated/markov_autoregression.html)

It seems I was more stupid than I thought. I'll update the codes soon.

@cc7768
Copy link
Member

cc7768 commented Oct 17, 2019

Excellent -- That sounds like a great place to find a test for the code :)

Many thanks for putting this together and for addressing the "nitpicks" that I have.

@Shunsuke-Hori
Copy link
Collaborator Author

I updated the code. Could anyone do another round of review?

Copy link
Member

@cc7768 cc7768 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this looks good now. Would be helpful if someone else took a review -- Maybe @sglyon can have a look and then merge when he's happy.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants